global_cycle  1.13.0
 All Data Structures Files Functions Variables Pages
cycle.f90
Go to the documentation of this file.
1 
4 
101  PROGRAM sfc_drv
102 
103  use mpi
104 
105  IMPLICIT NONE
106 !
107  CHARACTER(LEN=3) :: donst
108  INTEGER :: idim, jdim, lsoil, lugb, iy, im, id, ih, ialb
109  INTEGER :: isot, ivegsrc, lensfc, zsea1_mm, zsea2_mm, ierr
110  INTEGER :: nprocs, myrank, num_threads, num_parthds, max_tasks
111  REAL :: fh, deltsfc, zsea1, zsea2
112  LOGICAL :: use_ufo, do_nsst, do_lndinc, do_sfccycle, frac_grid
113 !
114  namelist/namcyc/ idim,jdim,lsoil,lugb,iy,im,id,ih,fh,&
115  deltsfc,ialb,use_ufo,donst, &
116  do_sfccycle,isot,ivegsrc,zsea1_mm, &
117  zsea2_mm, max_tasks, do_lndinc, frac_grid
118 !
119  DATA idim,jdim,lsoil/96,96,4/
120  DATA iy,im,id,ih,fh/1997,8,2,0,0./
121  DATA lugb/51/, deltsfc/0.0/, ialb/1/, max_tasks/99999/
122  DATA isot/1/, ivegsrc/2/, zsea1_mm/0/, zsea2_mm/0/
123 !
124  CALL mpi_init(ierr)
125  CALL mpi_comm_size(mpi_comm_world, nprocs, ierr)
126  CALL mpi_comm_rank(mpi_comm_world, myrank, ierr)
127 
128  if (myrank==0) call w3tagb('GLOBAL_CYCLE',2018,0179,0055,'NP20')
129 
130  num_threads = num_parthds()
131 
132  print*
133  print*,"STARTING CYCLE PROGRAM ON RANK ", myrank
134  print*,"RUNNING WITH ", nprocs, "TASKS"
135  print*,"AND WITH ", num_threads, " THREADS."
136 
137  use_ufo = .false.
138  donst = "NO"
139  do_lndinc = .false.
140  do_sfccycle = .true.
141  frac_grid = .false.
142 
143  print*
144  print*,"READ NAMCYC NAMELIST."
145 
146  CALL baopenr(36, "fort.36", ierr)
147  READ(36, nml=namcyc)
148 !IF (MYRANK==0) WRITE(6,NAMCYC)
149 
150  IF (max_tasks < 99999 .AND. myrank > (max_tasks - 1)) THEN
151  print*,"USER SPECIFIED MAX NUMBER OF TASKS: ", max_tasks
152  print*,"WILL NOT RUN CYCLE PROGRAM ON RANK: ", myrank
153  goto 333
154  ENDIF
155 
156  lensfc = idim*jdim ! TOTAL NUMBER OF POINTS FOR THE CUBED-SPHERE TILE
157 
158  zsea1 = float(zsea1_mm) / 1000.0 ! CONVERT FROM MM TO METERS
159  zsea2 = float(zsea2_mm) / 1000.0
160 
161  IF (donst == "YES") THEN
162  do_nsst=.true.
163  ELSE
164  do_nsst=.false.
165  ENDIF
166 
167  print*
168  IF (myrank==0) print*,"LUGB,IDIM,JDIM,ISOT,IVEGSRC,LSOIL,DELTSFC,IY,IM,ID,IH,FH: ", &
169  lugb,idim,jdim,isot,ivegsrc,lsoil,deltsfc,iy,im,id,ih,fh
170 
171  CALL sfcdrv(lugb,idim,jdim,lensfc,lsoil,deltsfc, &
172  iy,im,id,ih,fh,ialb, &
173  use_ufo,do_nsst,do_sfccycle,do_lndinc, &
174  frac_grid,zsea1,zsea2,isot,ivegsrc,myrank)
175 
176  print*
177  print*,'CYCLE PROGRAM COMPLETED NORMALLY ON RANK: ', myrank
178 
179  333 CONTINUE
180 
181  CALL mpi_barrier(mpi_comm_world, ierr)
182 
183  if (myrank==0) call w3tage('GLOBAL_CYCLE')
184 
185  CALL mpi_finalize(ierr)
186 
187  stop
188 
189  END PROGRAM sfc_drv
190 
299  SUBROUTINE sfcdrv(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, &
300  iy,im,id,ih,fh,ialb, &
301  use_ufo,do_nsst,do_sfccycle,do_lndinc,&
302  frac_grid,zsea1,zsea2,isot,ivegsrc,myrank)
303 !
304  USE read_write_data
305  use machine
306  USE mpi
307  USE land_increments, ONLY: add_increment_soil, &
308  add_increment_snow, &
309  calculate_landinc_mask, &
310  apply_land_da_adjustments_soil, &
311  apply_land_da_adjustments_snd, &
312  lsm_noah, lsm_noahmp
313 
314  IMPLICIT NONE
315 
316  INTEGER, INTENT(IN) :: idim, jdim, lensfc, lsoil, ialb
317  INTEGER, INTENT(IN) :: lugb, iy, im, id, ih
318  INTEGER, INTENT(IN) :: isot, ivegsrc, myrank
319 
320  LOGICAL, INTENT(IN) :: use_ufo, do_nsst,do_sfccycle
321  LOGICAL, INTENT(IN) :: do_lndinc, frac_grid
322 
323  REAL, INTENT(IN) :: fh, deltsfc, zsea1, zsea2
324 
325  INTEGER, PARAMETER :: nlunit=35
326  INTEGER, PARAMETER :: sz_nml=1
327 
328  CHARACTER(LEN=5) :: tile_num
329  CHARACTER(LEN=500) :: nst_file
330  CHARACTER(LEN=500) :: lnd_soi_file
331  CHARACTER(LEN=4) :: input_nml_file(sz_nml)
332 
333  INTEGER :: i, ierr
334  INTEGER :: i_index(lensfc), j_index(lensfc)
335  INTEGER :: idum(idim,jdim)
336  integer :: num_parthds, num_threads
337 
338  LOGICAL :: is_noahmp
339  INTEGER :: lsm
340 
341  real(kind=kind_io8) :: min_ice(lensfc)
342 
343  REAL :: slmask(lensfc), orog(lensfc)
344  REAL :: sihfcs(lensfc), sicfcs(lensfc)
345  REAL :: sitfcs(lensfc), tsffcs(lensfc)
346  REAL :: swefcs(lensfc), zorfcs(lensfc)
347  REAL :: albfcs(lensfc,4), tg3fcs(lensfc)
348  REAL :: cnpfcs(lensfc), smcfcs(lensfc,lsoil)
349  REAL :: stcfcs(lensfc,lsoil), slifcs(lensfc)
350  REAL :: aisfcs(lensfc), f10m(lensfc)
351  REAL :: vegfcs(lensfc), vetfcs(lensfc)
352  REAL :: sotfcs(lensfc), alffcs(lensfc,2)
353  REAL :: cvfcs(lensfc), cvtfcs(lensfc)
354  REAL :: cvbfcs(lensfc), tprcp(lensfc)
355  REAL :: srflag(lensfc), sndfcs(lensfc)
356  REAL :: slcfcs(lensfc,lsoil), vmxfcs(lensfc)
357  REAL :: vmnfcs(lensfc), t2m(lensfc)
358  REAL :: q2m(lensfc), slpfcs(lensfc)
359  REAL :: absfcs(lensfc), orog_uf(lensfc)
360  REAL :: ustar(lensfc)
361  REAL :: fmm(lensfc), fhh(lensfc)
362  REAL :: rla(lensfc), rlo(lensfc)
363  REAL(KIND=4) :: zsoil(lsoil)
364  REAL :: sig1t(lensfc)
367  REAL, ALLOCATABLE :: stc_bck(:,:), smc_bck(:,:), slc_bck(:,:)
368  REAL, ALLOCATABLE :: slifcs_fg(:), sicfcs_fg(:)
369  INTEGER, ALLOCATABLE :: landinc_mask_fg(:), landinc_mask(:)
370  REAL, ALLOCATABLE :: snd_bck(:), snd_inc(:), swe_bck(:)
371  REAL(KIND=KIND_IO8), ALLOCATABLE :: slmaskl(:), slmaskw(:), landfrac(:)
372 
373  TYPE(nsst_data) :: nsst
374  real, dimension(idim,jdim) :: tf_clm,tf_trd,sal_clm
375  real, dimension(lensfc) :: tf_clm_tile,tf_trd_tile,sal_clm_tile
376  INTEGER :: veg_type_landice
377  INTEGER, DIMENSION(LENSFC) :: stc_updated, slc_updated
378 
379  LOGICAL :: file_exists, do_soi_inc, do_sno_inc
380  CHARACTER(LEN=3) :: rankch
381 
382 !--------------------------------------------------------------------------------
383 ! NST_FILE is the path/name of the gaussian GSI file which contains NSST
384 ! increments.
385 !--------------------------------------------------------------------------------
386 
387  namelist/namsfcd/ nst_file, lnd_soi_file, do_sno_inc
388 
389  DATA nst_file/'NULL'/
390  DATA lnd_soi_file/'NULL'/
391 
392  do_sno_inc = .false.
393  do_soi_inc = .false.
394 
395  sig1t = 0.0 ! Not a dead start!
396 
397  input_nml_file = "NULL"
398 
399  CALL baopenr(37, "fort.37", ierr)
400  READ (37, nml=namsfcd)
401 
402  print*
403  print*,'IN ROUTINE SFCDRV,IDIM=',idim,'JDIM=',jdim,'FH=',fh
404 
405 !--------------------------------------------------------------------------------
406 ! READ THE OROGRAPHY AND GRID POINT LAT/LONS FOR THE CUBED-SPHERE TILE.
407 !--------------------------------------------------------------------------------
408 
409  ALLOCATE(landfrac(lensfc))
410  IF(frac_grid) THEN
411  print*,'- RUNNING WITH FRACTIONAL GRID.'
412  CALL read_lat_lon_orog(rla,rlo,orog,orog_uf,tile_num,idim,jdim,lensfc,landfrac=landfrac)
413  ELSE
414  CALL read_lat_lon_orog(rla,rlo,orog,orog_uf,tile_num,idim,jdim,lensfc)
415  landfrac=-999.9
416  ENDIF
417 
418  DO i = 1, idim
419  idum(i,:) = i
420  ENDDO
421 
422  i_index = reshape(idum, (/lensfc/))
423 
424  DO i = 1, jdim
425  idum(:,i) = i
426  ENDDO
427 
428  j_index = reshape(idum, (/lensfc/) )
429 
430  IF (do_nsst) THEN
431  print*
432  print*,"WILL PROCESS NSST RECORDS."
433  ALLOCATE(nsst%C_0(lensfc))
434  ALLOCATE(nsst%C_D(lensfc))
435  ALLOCATE(nsst%D_CONV(lensfc))
436  ALLOCATE(nsst%DT_COOL(lensfc))
437  ALLOCATE(nsst%IFD(lensfc))
438  ALLOCATE(nsst%QRAIN(lensfc))
439  ALLOCATE(nsst%TREF(lensfc))
440  ALLOCATE(nsst%TFINC(lensfc))
441  ALLOCATE(nsst%W_0(lensfc))
442  ALLOCATE(nsst%W_D(lensfc))
443  ALLOCATE(nsst%XS(lensfc))
444  ALLOCATE(nsst%XT(lensfc))
445  ALLOCATE(nsst%XTTS(lensfc))
446  ALLOCATE(nsst%XU(lensfc))
447  ALLOCATE(nsst%XV(lensfc))
448  ALLOCATE(nsst%XZ(lensfc))
449  ALLOCATE(nsst%XZTS(lensfc))
450  ALLOCATE(nsst%Z_C(lensfc))
451  ALLOCATE(nsst%ZM(lensfc))
452  ALLOCATE(slifcs_fg(lensfc))
453  ALLOCATE(sicfcs_fg(lensfc))
454  ENDIF
455 
456 IF (do_lndinc) THEN
457  ! identify variables to be updated, and allocate arrays.
458  IF (trim(lnd_soi_file) .NE. "NULL") THEN
459  do_soi_inc = .true.
460  print*
461  print*," APPLYING SOIL INCREMENTS FROM THE GSI"
462  ALLOCATE(stc_bck(lensfc, lsoil), smc_bck(lensfc, lsoil), slc_bck(lensfc,lsoil))
463  ALLOCATE(landinc_mask_fg(lensfc))
464  ENDIF
465  ! FOR NOW, CODE SO CAN DO BOTH, BUT MIGHT NEED TO THINK ABOUT THIS SOME MORE.
466  IF (do_sno_inc) THEN
467  ! ideally, would check here that sfcsub snow DA update is not also requested
468  ! but latter is controlled by fnsol, which is read in within that routine.
469  ! should be done at script level.
470  print*
471  print*," APPLYING SNOW INCREMENTS FROM JEDI"
472  ALLOCATE(snd_bck(lensfc), snd_inc(lensfc), swe_bck(lensfc))
473  ENDIF
474  ! set-up land mask info
475  ALLOCATE(landinc_mask(lensfc))
476  if (ivegsrc == 2) then ! sib
477  veg_type_landice=13
478  else
479  veg_type_landice=15
480  endif
481 ENDIF
482 
483 !--------------------------------------------------------------------------------
484 ! READ THE INPUT SURFACE DATA ON THE CUBED-SPHERE TILE.
485 !--------------------------------------------------------------------------------
486 
487  CALL read_data(lsoil,lensfc,do_nsst,.false.,is_noahmp=is_noahmp, &
488  tsffcs=tsffcs,smcfcs=smcfcs, &
489  swefcs=swefcs,stcfcs=stcfcs,tg3fcs=tg3fcs,zorfcs=zorfcs, &
490  cvfcs=cvfcs, cvbfcs=cvbfcs,cvtfcs=cvtfcs,albfcs=albfcs, &
491  vegfcs=vegfcs,slifcs=slifcs,cnpfcs=cnpfcs,f10m=f10m , &
492  vetfcs=vetfcs,sotfcs=sotfcs,alffcs=alffcs,ustar=ustar , &
493  fmm=fmm ,fhh=fhh ,sihfcs=sihfcs,sicfcs=sicfcs, &
494  sitfcs=sitfcs,tprcp=tprcp ,srflag=srflag,sndfcs=sndfcs, &
495  vmnfcs=vmnfcs,vmxfcs=vmxfcs,slcfcs=slcfcs,slpfcs=slpfcs, &
496  absfcs=absfcs,t2m=t2m ,q2m=q2m ,slmask=slmask, &
497  zsoil=zsoil, nsst=nsst)
498 
499  IF (frac_grid .AND. .NOT. is_noahmp) THEN
500  print *, 'FATAL ERROR: NOAH lsm update does not work with fractional grids.'
501  call mpi_abort(mpi_comm_world, 18, ierr)
502  ENDIF
503 
504  IF (is_noahmp .AND. do_sno_inc) THEN
505  print *, 'FATAL ERROR: Snow increment update does not work with NOAH_MP.'
506  call mpi_abort(mpi_comm_world, 29, ierr)
507  ENDIF
508 
509  IF (is_noahmp) THEN
510  lsm=lsm_noahmp
511  ELSE
512  lsm=lsm_noah
513  ENDIF
514 
515  IF (use_ufo) THEN
516  print*
517  print*,'USE UNFILTERED OROGRAPHY.'
518  ELSE
519  orog_uf = 0.0
520  ENDIF
521 
522  DO i=1,lensfc
523  aisfcs(i) = 0.
524  IF(nint(slifcs(i)).EQ.2) aisfcs(i) = 1.
525  ENDDO
526 
527  IF (do_nsst) THEN
528  sicfcs_fg=sicfcs
529  IF (.NOT. do_sfccycle ) THEN
530  print*
531  print*,"FIRST GUESS MASK ADJUSTED BY IFD RECORD"
532  slifcs_fg = slifcs
533  WHERE(nint(nsst%IFD) == 3) slifcs_fg = 2.0
534  ELSE
535  print*
536  print*,"SAVE FIRST GUESS MASK"
537  slifcs_fg = slifcs
538  ENDIF
539  ENDIF
540 
541  ! CALCULATE MASK FOR LAND INCREMENTS
542  IF (do_lndinc) &
543  CALL calculate_landinc_mask(smcfcs(:,1),swefcs, vetfcs, &
544  lensfc,veg_type_landice, landinc_mask)
545 
546 !--------------------------------------------------------------------------------
547 ! UPDATE SURFACE FIELDS.
548 !
549 ! FIRST, SET WATER AND LAND MASKS - SLMASKW/SLMASKL. FOR UNCOUPLED
550 ! (NON-FRACTIONAL) MODE, THESE ARE IDENTICAL TO THE CURRENT
551 ! MASK - '0' WATER; '1' LAND.
552 !--------------------------------------------------------------------------------
553 
554  IF (do_sfccycle) THEN
555  ALLOCATE(slmaskl(lensfc), slmaskw(lensfc))
556 
557  set_mask : IF (frac_grid) THEN
558 
559  DO i=1,lensfc
560  IF(landfrac(i) > 0.0_kind_io8) THEN
561  slmaskl(i) = ceiling(landfrac(i)-1.0e-6_kind_io8)
562  slmaskw(i) = floor(landfrac(i)+1.0e-6_kind_io8)
563  ELSE
564  IF(nint(slmask(i)) == 1) THEN ! If landfrac is zero, this should not happen.
565  ! So, stop processing.
566  print*, 'FATAL ERROR: LAND FRAC AND SLMASK MISMATCH.'
567  CALL mpi_abort(mpi_comm_world, 27, ierr)
568  ELSE
569  slmaskl(i) = 0.0_kind_io8
570  slmaskw(i) = 0.0_kind_io8
571  ENDIF
572  ENDIF
573 
574  ENDDO
575 
576  ELSE
577 
578 ! For running uncoupled (non-fractional grid).
579 
580  DO i=1,lensfc
581  IF(nint(slmask(i)) == 1) THEN
582  slmaskl(i) = 1.0_kind_io8
583  slmaskw(i) = 1.0_kind_io8
584  ELSE
585  slmaskl(i) = 0.0_kind_io8
586  slmaskw(i) = 0.0_kind_io8
587  ENDIF
588  ENDDO
589 
590  ENDIF set_mask
591 
592  DO i=1,lensfc
593  if(nint(slmask(i)) == 0) then
594  min_ice(i) = 0.15_kind_io8
595  else
596  min_ice(i) = 0.0_kind_io8
597  endif
598  ENDDO
599 
600  num_threads = num_parthds()
601  print*
602  print*,"CALL SFCCYCLE TO UPDATE SURFACE FIELDS."
603  CALL sfccycle(lugb,lensfc,lsoil,sig1t,deltsfc, &
604  iy,im,id,ih,fh,rla,rlo, &
605  slmaskl,slmaskw, orog, orog_uf, use_ufo, do_nsst, &
606  sihfcs,sicfcs,sitfcs,sndfcs,slcfcs, &
607  vmnfcs,vmxfcs,slpfcs,absfcs, &
608  tsffcs,swefcs,zorfcs,albfcs,tg3fcs, &
609  cnpfcs,smcfcs,stcfcs,slifcs,aisfcs, &
610  vegfcs,vetfcs,sotfcs,alffcs, &
611  cvfcs,cvbfcs,cvtfcs,myrank,num_threads, nlunit, &
612  sz_nml, input_nml_file, &
613  min_ice, &
614  ialb,isot,ivegsrc,tile_num,i_index,j_index)
615 
616  DEALLOCATE(slmaskl, slmaskw)
617  ENDIF
618 
619 !--------------------------------------------------------------------------------
620 ! IF RUNNING WITH NSST, READ IN GSI FILE WITH THE UPDATED INCREMENTS (ON THE
621 ! GAUSSIAN GRID), INTERPOLATE INCREMENTS TO THE CUBED-SPHERE TILE, AND PERFORM
622 ! REQUIRED ADJUSTMENTS AND QC.
623 !--------------------------------------------------------------------------------
624 
625  IF (do_nsst) THEN
626  IF (nst_file == "NULL") THEN
627  print*
628  print*,"NO GSI FILE. ADJUST IFD FOR FORMER ICE POINTS."
629  DO i = 1, lensfc
630  IF (sicfcs_fg(i) > 0.0 .AND. sicfcs(i) == 0) THEN
631  nsst%IFD(i) = 3.0
632  ENDIF
633  ENDDO
634  nsst%TFINC = 0.0
635  ELSE
636  print*
637  print*,"ADJUST TREF FROM GSI INCREMENT"
638 !
639 ! Get tf climatology at the time
640 !
641  call get_tf_clm(rla,rlo,jdim,idim,iy,im,id,ih,tf_clm,tf_trd)
642  tf_clm_tile(:) = reshape(tf_clm, (/lensfc/) )
643  tf_trd_tile(:) = reshape(tf_trd, (/lensfc/) )
644 !
645 ! Get salinity climatology at the time
646 !
647  call get_sal_clm(rla,rlo,jdim,idim,iy,im,id,ih,sal_clm)
648  sal_clm_tile(:) = reshape(sal_clm, (/lensfc/) )
649 !
650 ! read tf analysis increment generated by GSI
651 !
652  CALL read_gsi_data(nst_file, 'NST')
653 !
654 ! update foundation & surface temperature for NSST
655 !
656  CALL adjust_nsst(rla,rlo,slifcs,slifcs_fg,tsffcs,sitfcs,sicfcs,sicfcs_fg,&
657  stcfcs,nsst,lensfc,lsoil,idim,jdim,zsea1,zsea2, &
658  tf_clm_tile,tf_trd_tile,sal_clm_tile,landfrac,frac_grid)
659  ENDIF
660  ENDIF
661 
662 !--------------------------------------------------------------------------------
663 ! READ IN AND APPLY LAND INCREMENTS FROM THE GSI
664 !--------------------------------------------------------------------------------
665 
666  IF (do_lndinc) THEN
667 
668  ! SNOW INCREMENTS
669  ! do snow first, as temperature updates will use snow analaysis
670  IF (do_sno_inc) THEN
671  ! updates are made to snow depth only over land (and not-land ice).
672  ! SWE is then updated from the snow depth analysis, using the model
673  ! forecast density
674 
675  !--------------------------------------------------------------------------------
676  ! read increments in
677  !--------------------------------------------------------------------------------
678 
679  ! Only coded for DA on native model grid (would always be the case for cycling DA)
680  CALL read_data(lsoil,lensfc,.false.,.true.,sndfcs=snd_inc)
681 
682  !--------------------------------------------------------------------------------
683  ! add increments to state vars
684  !--------------------------------------------------------------------------------
685 
686  ! store background states
687  snd_bck = sndfcs
688  swe_bck = swefcs
689 
690  CALL add_increment_snow(snd_inc,landinc_mask,lensfc,sndfcs)
691 
692  !--------------------------------------------------------------------------------
693  ! make any necessary adjustments to dependent variables
694  !--------------------------------------------------------------------------------
695 
696  CALL apply_land_da_adjustments_snd(lsm, lensfc, landinc_mask, swe_bck, snd_bck, &
697  sndfcs, swefcs)
698 
699  ENDIF
700 
701  ! SOIL INCREMENTS
702  IF (do_soi_inc) THEN
703 
704  !--------------------------------------------------------------------------------
705  ! re-calculate soilsnow mask if snow has been updated.
706  !--------------------------------------------------------------------------------
707 
708  landinc_mask_fg = landinc_mask
709 
710  IF (do_sfccycle .OR. do_sno_inc) THEN
711  CALL calculate_landinc_mask(smcfcs(:,1),swefcs, vetfcs, lensfc, &
712  veg_type_landice, landinc_mask )
713  ENDIF
714 
715  !--------------------------------------------------------------------------------
716  ! read increments in
717  !--------------------------------------------------------------------------------
718 
719  WRITE(rankch, '(I3.3)') (myrank+1)
720 
721  lnd_soi_file = trim(lnd_soi_file) // "." // rankch
722 
723  INQUIRE(file=trim(lnd_soi_file), exist=file_exists)
724  IF (.not. file_exists) then
725  print *, 'FATAL ERROR: land increment update requested, but file does not exist: ', &
726  trim(lnd_soi_file)
727  call mpi_abort(mpi_comm_world, 10, ierr)
728  ENDIF
729 
730  CALL read_gsi_data(lnd_soi_file, 'LND', lsoil=lsoil)
731 
732  !--------------------------------------------------------------------------------
733  ! add increments to state vars
734  !--------------------------------------------------------------------------------
735  ! when applying increments, will often need to adjust other land states in response
736  ! to the changes made. Need to store bacground, apply the increments, then make
737  ! secondart adjustments. When updating more than one state, be careful about the
738  ! order if increments and secondary adjustments.
739 
740  ! store background states
741  stc_bck = stcfcs
742  smc_bck = smcfcs
743  slc_bck = slcfcs
744 
745  ! below updates [STC/SMC/STC]FCS to hold the analysis
746  CALL add_increment_soil(rla,rlo,stcfcs,smcfcs,slcfcs,stc_updated, slc_updated, &
747  landinc_mask,landinc_mask_fg,lensfc,lsoil,idim,jdim,lsm,myrank)
748 
749  !--------------------------------------------------------------------------------
750  ! make any necessary adjustments to dependent variables
751  !--------------------------------------------------------------------------------
752 
753 
754  CALL apply_land_da_adjustments_soil(lsm, isot, ivegsrc,lensfc, lsoil, &
755  sotfcs, landinc_mask_fg, stc_bck, stcfcs, smcfcs, slcfcs, stc_updated, &
756  slc_updated,zsoil)
757 
758  ENDIF ! soil increments
759 
760 !--------------------------------------------------------------------------------
761 ! clean up
762 !--------------------------------------------------------------------------------
763 
764  ! to do - save and write out STC_INC? (soil temperature increments)
765  IF(ALLOCATED(landinc_mask_fg)) DEALLOCATE(landinc_mask_fg)
766  IF(ALLOCATED(landinc_mask)) DEALLOCATE(landinc_mask)
767  IF(ALLOCATED(stc_bck)) DEALLOCATE(stc_bck)
768  IF(ALLOCATED(smc_bck)) DEALLOCATE(smc_bck)
769  IF(ALLOCATED(slc_bck)) DEALLOCATE(slc_bck)
770  IF(ALLOCATED(snd_bck)) DEALLOCATE(snd_bck)
771  IF(ALLOCATED(swe_bck)) DEALLOCATE(swe_bck)
772  IF(ALLOCATED(snd_inc)) DEALLOCATE(snd_inc)
773 
774  ENDIF
775 !--------------------------------------------------------------------------------
776 ! WRITE OUT UPDATED SURFACE DATA ON THE CUBED-SPHERE TILE.
777 !--------------------------------------------------------------------------------
778 
779  IF (lsm==lsm_noahmp) THEN
780 
781  CALL write_data(lensfc,idim,jdim,lsoil,do_nsst,nsst,vegfcs=vegfcs, &
782  slcfcs=slcfcs,smcfcs=smcfcs,stcfcs=stcfcs,&
783  sicfcs=sicfcs,sihfcs=sihfcs)
784 
785  ELSEIF (lsm==lsm_noah) THEN
786 
787  CALL write_data(lensfc,idim,jdim,lsoil, &
788  do_nsst,nsst,slifcs=slifcs,tsffcs=tsffcs,vegfcs=vegfcs, &
789  swefcs=swefcs,tg3fcs=tg3fcs,zorfcs=zorfcs, &
790  albfcs=albfcs,alffcs=alffcs,cnpfcs=cnpfcs, &
791  f10m=f10m,t2m=t2m,q2m=q2m,vetfcs=vetfcs, &
792  sotfcs=sotfcs,ustar=ustar,fmm=fmm,fhh=fhh, &
793  sicfcs=sicfcs,sihfcs=sihfcs,sitfcs=sitfcs,tprcp=tprcp, &
794  srflag=srflag,swdfcs=sndfcs,vmnfcs=vmnfcs, &
795  vmxfcs=vmxfcs,slpfcs=slpfcs,absfcs=absfcs, &
796  slcfcs=slcfcs,smcfcs=smcfcs,stcfcs=stcfcs)
797 
798  ENDIF
799 
800  IF (do_nsst) THEN
801  DEALLOCATE(nsst%C_0)
802  DEALLOCATE(nsst%C_D)
803  DEALLOCATE(nsst%D_CONV)
804  DEALLOCATE(nsst%DT_COOL)
805  DEALLOCATE(nsst%IFD)
806  DEALLOCATE(nsst%QRAIN)
807  DEALLOCATE(nsst%TREF)
808  DEALLOCATE(nsst%TFINC)
809  DEALLOCATE(nsst%W_0)
810  DEALLOCATE(nsst%W_D)
811  DEALLOCATE(nsst%XS)
812  DEALLOCATE(nsst%XT)
813  DEALLOCATE(nsst%XTTS)
814  DEALLOCATE(nsst%XU)
815  DEALLOCATE(nsst%XV)
816  DEALLOCATE(nsst%XZ)
817  DEALLOCATE(nsst%XZTS)
818  DEALLOCATE(nsst%Z_C)
819  DEALLOCATE(nsst%ZM)
820  DEALLOCATE(slifcs_fg)
821  DEALLOCATE(sicfcs_fg)
822  ENDIF
823 
824  RETURN
825 
826  END SUBROUTINE sfcdrv
827 
859  SUBROUTINE adjust_nsst(RLA,RLO,SLMSK_TILE,SLMSK_FG_TILE,SKINT_TILE,&
860  sicet_tile,sice_tile,sice_fg_tile,soilt_tile,nsst, &
861  lensfc,lsoil,idim,jdim,zsea1,zsea2, &
862  tf_clm_tile,tf_trd_tile,sal_clm_tile,landfrac, &
863  frac_grid)
864 
865  USE utils
866  USE gdswzd_mod
867  USE read_write_data, ONLY : idim_gaus, jdim_gaus, &
868  slmsk_gaus, dtref_gaus, &
869  nsst_data
870 
871  USE mpi
872 
873  IMPLICIT NONE
874 
875  INTEGER, INTENT(IN) :: lensfc, lsoil, idim, jdim
876 
877  LOGICAL, INTENT(IN) :: frac_grid
878 
879  REAL, INTENT(IN) :: slmsk_tile(lensfc), slmsk_fg_tile(lensfc), landfrac(lensfc)
880  real, intent(in) :: tf_clm_tile(lensfc),tf_trd_tile(lensfc),sal_clm_tile(lensfc)
881  REAL, INTENT(IN) :: zsea1, zsea2,sice_tile(lensfc),sice_fg_tile(lensfc)
882  REAL, INTENT(IN) :: rla(lensfc), rlo(lensfc)
883  REAL, INTENT(INOUT) :: skint_tile(lensfc)
884  REAL, INTENT(INOUT) :: sicet_tile(lensfc),soilt_tile(lensfc,lsoil)
885 
886  TYPE(nsst_data) :: nsst
887 
888  REAL, PARAMETER :: tmax=313.0,tzero=273.16
889 
890  INTEGER :: iopt, nret, kgds_gaus(200)
891  INTEGER :: igaus, jgaus, ij, ii, jj, iii, jjj, krad
892  INTEGER :: istart, iend, jstart, jend
893 !INTEGER :: MASK_TILE, MASK_FG_TILE
894  INTEGER,allocatable :: mask_tile(:),mask_fg_tile(:)
895  INTEGER :: itile, jtile
896  INTEGER :: max_search, j, ierr
897  INTEGER :: igausp1, jgausp1
898  integer :: nintp,nsearched,nice,nland
899  integer :: nfill,nfill_tice,nfill_clm
900  integer :: nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c
901 
902  INTEGER, ALLOCATABLE :: id1(:,:), id2(:,:), jdc(:,:)
903 
904  LOGICAL :: is_ice
905 
906  real :: tfreez
907  REAL :: tref_save,wsum,tf_ice,tf_thaw
908  REAL :: fill, dtzm, gaus_res_km, dtref
909  REAL, ALLOCATABLE :: xpts(:), ypts(:), lats(:), lons(:)
910  REAL, ALLOCATABLE :: dum2d(:,:), lats_rad(:), lons_rad(:)
911  REAL, ALLOCATABLE :: agrid(:,:,:), s2c(:,:,:)
912 
913  kgds_gaus = 0
914  kgds_gaus(1) = 4 ! OCT 6 - TYPE OF GRID (GAUSSIAN)
915  kgds_gaus(2) = idim_gaus ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
916  kgds_gaus(3) = jdim_gaus
917  kgds_gaus(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN
918  kgds_gaus(5) = 0 ! OCT 14-16 - LON OF ORIGIN
919  kgds_gaus(6) = 128 ! OCT 17 - RESOLUTION FLAG
920  kgds_gaus(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT
921  kgds_gaus(8) = nint(-360000./float(idim_gaus)) ! OCT 21-23 - LON OF EXTREME POINT
922  kgds_gaus(9) = nint((360.0 / float(idim_gaus))*1000.0)
923  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
924  kgds_gaus(10) = jdim_gaus/2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
925  kgds_gaus(12) = 255 ! OCT 29 - RESERVED
926  kgds_gaus(20) = 255 ! OCT 5 - NOT USED, SET TO 255
927 
928  print*
929  print*,'ADJUST NSST USING GSI INCREMENTS ON GAUSSIAN GRID'
930 
931 !----------------------------------------------------------------------
932 ! CALL GDSWZD TO COMPUTE THE LAT/LON OF EACH GSI GAUSSIAN GRID POINT.
933 !----------------------------------------------------------------------
934 
935  iopt = 0
936  fill = -9999.
937  ALLOCATE(xpts(idim_gaus*jdim_gaus))
938  ALLOCATE(ypts(idim_gaus*jdim_gaus))
939  ALLOCATE(lats(idim_gaus*jdim_gaus))
940  ALLOCATE(lons(idim_gaus*jdim_gaus))
941  xpts = fill
942  ypts = fill
943  lats = fill
944  lons = fill
945 
946  CALL gdswzd(kgds_gaus,iopt,(idim_gaus*jdim_gaus),fill,xpts,ypts,lons,lats,nret)
947 
948  IF (nret /= (idim_gaus*jdim_gaus)) THEN
949  print*,'FATAL ERROR: PROBLEM IN GDSWZD. STOP.'
950  CALL mpi_abort(mpi_comm_world, 12, ierr)
951  ENDIF
952 
953  DEALLOCATE (xpts, ypts)
954 
955  ALLOCATE(dum2d(idim_gaus,jdim_gaus))
956  dum2d = reshape(lats, (/idim_gaus,jdim_gaus/) )
957  DEALLOCATE(lats)
958 
959  ALLOCATE(lats_rad(jdim_gaus))
960  DO j = 1, jdim_gaus
961  lats_rad(j) = dum2d(1,jdim_gaus-j+1) * 3.1415926 / 180.0
962  ENDDO
963 
964  dum2d = reshape(lons, (/idim_gaus,jdim_gaus/) )
965  DEALLOCATE(lons)
966  ALLOCATE(lons_rad(idim_gaus))
967  lons_rad = dum2d(:,1) * 3.1415926 / 180.0
968 
969  DEALLOCATE(dum2d)
970 
971  ALLOCATE(agrid(idim,jdim,2))
972  agrid(:,:,1) = reshape(rlo, (/idim,jdim/) )
973  agrid(:,:,2) = reshape(rla, (/idim,jdim/) )
974  agrid = agrid * 3.1415926 / 180.0
975 
976  ALLOCATE(id1(idim,jdim))
977  ALLOCATE(id2(idim,jdim))
978  ALLOCATE(jdc(idim,jdim))
979  ALLOCATE(s2c(idim,jdim,4))
980 
981 !----------------------------------------------------------------------
982 ! COMPUTE BILINEAR WEIGHTS FOR EACH MODEL POINT FROM THE NEAREST
983 ! FOUR GSI/GAUSSIAN POINTS. DOES NOT ACCOUNT FOR MASK. THAT
984 ! HAPPENS LATER.
985 !----------------------------------------------------------------------
986 
987  CALL remap_coef( 1, idim, 1, jdim, idim_gaus, jdim_gaus, &
988  lons_rad, lats_rad, id1, id2, jdc, s2c, agrid )
989 
990  DEALLOCATE(lons_rad, lats_rad, agrid)
991 
992 !----------------------------------------------------------------------
993 ! THE MAXIMUM DISTANCE TO SEARCH IS 500 KM. HOW MANY GAUSSIAN
994 ! GRID LENGTHS IS THAT?
995 !----------------------------------------------------------------------
996 
997  gaus_res_km = 360.0 / idim_gaus * 111.0
998  max_search = ceiling(500.0/gaus_res_km)
999 
1000  print*
1001  print*,'MAXIMUM SEARCH IS ',max_search, ' GAUSSIAN POINTS.'
1002  print*
1003 
1004 !
1005 ! Initialize variables for counts statitics to be zeros
1006 !
1007  nintp = 0
1008  nset_thaw = 0
1009  nset_thaw_s = 0
1010  nset_thaw_i = 0
1011  nset_thaw_c = 0
1012  nsearched = 0
1013  nfill = 0
1014  nfill_tice = 0
1015  nfill_clm = 0
1016  nice = 0
1017  nland = 0
1018 !----------------------------------------------------------------------
1019 ! TREF INCREMENT WILL BE OUTPUT. INITIALIZE TO ZERO.
1020 !----------------------------------------------------------------------
1021 
1022  nsst%TFINC = 0.0
1023 
1024  allocate(mask_tile(lensfc))
1025  allocate(mask_fg_tile(lensfc))
1026 
1027  IF(.NOT. frac_grid) THEN
1028  mask_tile = nint(slmsk_tile)
1029  mask_fg_tile = nint(slmsk_fg_tile)
1030  ELSE
1031  mask_tile=0
1032  WHERE(sice_tile > 0.0) mask_tile=2
1033  WHERE(landfrac == 1.0) mask_tile=1
1034  mask_fg_tile=0
1035  WHERE(sice_fg_tile > 0.0) mask_fg_tile=2
1036  WHERE(landfrac == 1.0) mask_fg_tile=1
1037  ENDIF
1038 
1039  ij_loop : DO ij = 1, lensfc
1040 
1041 !
1042 ! when sea ice exists, get salinity dependent water temperature
1043 !
1044  tf_ice = tfreez(sal_clm_tile(ij)) + tzero
1045 !----------------------------------------------------------------------
1046 ! SKIP LAND POINTS. NSST NOT APPLIED AT LAND.
1047 !----------------------------------------------------------------------
1048 
1049  IF (mask_tile(ij) == 1) THEN
1050  nland = nland + 1
1051  cycle ij_loop
1052  ENDIF
1053 
1054 !
1055 ! these are ice points. set tref to tf_ice and update tmpsfc.
1056 !
1057  if (mask_tile(ij) == 2) then
1058  nsst%tref(ij)=tf_ice ! water part tmp set
1059  skint_tile(ij)=(1.0-sice_tile(ij))*nsst%tref(ij)+sice_tile(ij)*sicet_tile(ij)
1060  nice = nice + 1
1061  cycle ij_loop
1062  endif
1063 
1064 !
1065 ! Get i,j index on array of (idim,jdim) from known ij
1066 !
1067  jtile = (ij-1) / idim + 1
1068  itile = mod(ij,idim)
1069  IF (itile==0) itile = idim
1070 
1071 !----------------------------------------------------------------------
1072 ! IF THE MODEL POINT WAS ICE COVERED, BUT IS NOW OPEN WATER, SET
1073 ! TREF TO searched adjascent open water onea, if failed the search, set to
1074 ! weighted average of tf_ice and tf_clm. For NSST vars, set xz TO '30' AND ALL OTHER FIELDS TO ZERO.
1075 !----------------------------------------------------------------------
1076 
1077  IF (mask_fg_tile(ij) == 2 .AND. mask_tile(ij) == 0) THEN
1078 !
1079 ! set background for the thaw (just melted water) situation
1080 !
1081  call tf_thaw_set(nsst%tref,mask_fg_tile,itile,jtile,tf_ice,tf_clm_tile(ij),tf_thaw,idim,jdim, &
1082  nset_thaw_s,nset_thaw_i,nset_thaw_c)
1083  call nsst_water_reset(nsst,ij,tf_thaw)
1084  nset_thaw = nset_thaw + 1
1085  ENDIF
1086 
1087 !----------------------------------------------------------------------
1088 ! THESE ARE POINTS THAT ARE OPEN WATER AND WERE OPEN WATER PRIOR
1089 ! TO ANY ICE UPDATE BY SFCCYCLE. UPDATE TREF AND SKIN TEMP.
1090 ! AT OPEN WATER POINTS, THE SEA ICE TEMPERATURE (SICET_TILE) AND
1091 ! SOIL COLUMN TEMPERATURE (SOILT_TILE) ARE SET TO THE SKIN TEMP.
1092 ! IT IS SIMPLY A FILLER VALUE. THESE FIELDS ARE NOT USED AT
1093 ! OPEN WATER POINTS.
1094 !----------------------------------------------------------------------
1095 !----------------------------------------------------------------------
1096 ! SEE IF ANY OF THE NEAREST GSI POINTS MASK AREA OPEN WATER.
1097 ! IF SO, APPLY NSST INCREMENT USING BILINEAR INTERPOLATION.
1098 !----------------------------------------------------------------------
1099 
1100  igaus = id1(itile,jtile)
1101  jgaus = jdc(itile,jtile)
1102  igausp1 = id2(itile,jtile)
1103  jgausp1 = jdc(itile,jtile)+1
1104 
1105  IF (slmsk_gaus(igaus,jgaus) == 0 .OR. &
1106  slmsk_gaus(igausp1,jgaus) == 0 .OR. &
1107  slmsk_gaus(igausp1,jgausp1) == 0 .OR. &
1108  slmsk_gaus(igaus,jgausp1) == 0) THEN
1109 
1110  dtref = 0.0
1111  wsum = 0.0
1112 
1113  IF (slmsk_gaus(igaus,jgaus) == 0) THEN
1114  dtref = dtref + (s2c(itile,jtile,1) * dtref_gaus(igaus,jgaus))
1115  wsum = wsum + s2c(itile,jtile,1)
1116  ENDIF
1117 
1118  IF (slmsk_gaus(igausp1,jgaus) == 0) THEN
1119  dtref = dtref + (s2c(itile,jtile,2) * dtref_gaus(igausp1,jgaus))
1120  wsum = wsum + s2c(itile,jtile,2)
1121  ENDIF
1122 
1123  IF (slmsk_gaus(igausp1,jgausp1) == 0) THEN
1124  dtref = dtref + (s2c(itile,jtile,3) * dtref_gaus(igausp1,jgausp1))
1125  wsum = wsum + s2c(itile,jtile,3)
1126  ENDIF
1127 
1128  IF (slmsk_gaus(igaus,jgausp1) == 0) THEN
1129  dtref = dtref + (s2c(itile,jtile,4) * dtref_gaus(igaus,jgausp1))
1130  wsum = wsum + s2c(itile,jtile,4)
1131  ENDIF
1132 
1133  nintp = nintp + 1
1134  dtref = dtref / wsum
1135 
1136  tref_save = nsst%TREF(ij)
1137  nsst%TREF(ij) = nsst%TREF(ij) + dtref
1138  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
1139  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
1140  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
1141 
1142  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
1143  nsst%Z_C(ij),zsea1,zsea2,dtzm)
1144 
1145  skint_tile(ij) = nsst%TREF(ij) + dtzm
1146  skint_tile(ij) = max(skint_tile(ij), tf_ice)
1147  skint_tile(ij) = min(skint_tile(ij), tmax)
1148 
1149  sicet_tile(ij) = skint_tile(ij)
1150 ! Under fractional grids, soilt is used at points with at
1151 ! least some land.
1152  IF(.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)
1153 
1154 !----------------------------------------------------------------------
1155 ! NO NEARBY GSI/GAUSSIAN OPEN WATER POINTS. PERFORM A SPIRAL SEARCH TO
1156 ! FIND NEAREST NON-LAND POINT ON GSI/GAUSSIAN GRID.
1157 !----------------------------------------------------------------------
1158 
1159  ELSE
1160 
1161  is_ice = .false.
1162 
1163  DO krad = 1, max_search
1164 
1165  istart = igaus - krad
1166  iend = igaus + krad
1167  jstart = jgaus - krad
1168  jend = jgaus + krad
1169 
1170  DO jj = jstart, jend
1171  DO ii = istart, iend
1172 
1173  IF((jj == jstart) .OR. (jj == jend) .OR. &
1174  (ii == istart) .OR. (ii == iend)) THEN
1175 
1176  IF ((jj >= 1) .AND. (jj <= jdim_gaus)) THEN
1177 
1178  jjj = jj
1179  IF (ii <= 0) THEN
1180  iii = idim_gaus + ii
1181  ELSE IF (ii >= (idim_gaus+1)) THEN
1182  iii = ii - idim_gaus
1183  ELSE
1184  iii = ii
1185  END IF
1186 
1187 !----------------------------------------------------------------------
1188 ! SEE IF NEARBY POINTS ARE SEA ICE. IF THEY ARE, AND THE SEARCH FOR
1189 ! A GAUSSIAN GRID OPEN WATER POINT FAILS, THEN TREF WILL BE SET TO
1190 ! FREEZING BELOW.
1191 !----------------------------------------------------------------------
1192 
1193  IF (krad <= 2 .AND. slmsk_gaus(iii,jjj) == 2) is_ice = .true.
1194 
1195  IF (slmsk_gaus(iii,jjj) == 0) THEN
1196 
1197 ! PRINT*,'MISMATCH AT TILE POINT ',ITILE,JTILE
1198 ! PRINT*,'UPDATE TREF USING GSI INCREMENT AT ',III,JJJ,DTREF_GAUS(III,JJJ)
1199  nsearched = nsearched + 1
1200 
1201  tref_save = nsst%TREF(ij)
1202  nsst%TREF(ij ) = nsst%TREF(ij) + dtref_gaus(iii,jjj)
1203  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
1204  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
1205  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
1206 
1207  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
1208  nsst%Z_C(ij),zsea1,zsea2,dtzm)
1209 
1210  skint_tile(ij) = nsst%TREF(ij) + dtzm
1211  skint_tile(ij) = max(skint_tile(ij), tf_ice)
1212  skint_tile(ij) = min(skint_tile(ij), tmax)
1213 
1214  sicet_tile(ij) = skint_tile(ij)
1215  IF(.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)
1216  cycle ij_loop
1217 
1218  ENDIF ! GSI/Gaussian mask is open water
1219 
1220  ENDIF
1221 
1222  ENDIF
1223 
1224  ENDDO
1225  ENDDO
1226 
1227  ENDDO ! KRAD LOOP
1228 
1229 !----------------------------------------------------------------------
1230 ! THE SEARCH FAILED. IF THERE IS NEARBY ICE, SET TREF TO FREEZING.
1231 ! ELSE UPDATE TREF BASED ON THE ANNUAL SST CYCLE.
1232 !----------------------------------------------------------------------
1233 
1234 ! PRINT*,'WARNING !!!!!! SEARCH FAILED AT TILE POINT ',ITILE,JTILE
1235 
1236  nfill = nfill + 1
1237  IF (is_ice) THEN
1238  nsst%TREF(ij) = tf_ice
1239 ! PRINT*,"NEARBY ICE. SET TREF TO FREEZING"
1240  nfill_tice = nfill_tice + 1
1241  ELSE
1242  tref_save = nsst%TREF(ij)
1243  nsst%TREF(ij) = nsst%TREF(ij) + tf_trd_tile(ij)
1244  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
1245  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
1246  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
1247 ! PRINT*,'UPDATE TREF FROM SST CLIMO ',DTREF
1248  nfill_clm = nfill_clm + 1
1249  ENDIF
1250 
1251  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
1252  nsst%Z_C(ij),zsea1,zsea2,dtzm)
1253 
1254  skint_tile(ij) = nsst%TREF(ij) + dtzm
1255  skint_tile(ij) = max(skint_tile(ij), tf_ice)
1256  skint_tile(ij) = min(skint_tile(ij), tmax)
1257 
1258  sicet_tile(ij) = skint_tile(ij)
1259  IF (.NOT. frac_grid) soilt_tile(ij,:) = skint_tile(ij)
1260 
1261  ENDIF ! NEARBY GAUSSIAN POINTS ARE OPEN WATER?
1262 
1263  ENDDO ij_loop
1264 
1265  write(*,'(a)') 'statistics of grids number processed for tile : '
1266  write(*,'(a,I8)') ' nintp = ',nintp
1267  write(*,'(a,4I8)') 'nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c =',nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c
1268  write(*,'(a,I8)') ' nsearched = ',nsearched
1269  write(*,'(a,3I6)') ' nfill,nfill_tice,nfill_clm = ',nfill,nfill_tice,nfill_clm
1270  write(*,'(a,I8)') ' nice = ',nice
1271  write(*,'(a,I8)') ' nland = ',nland
1272 
1273  DEALLOCATE(id1, id2, jdc, s2c, mask_tile, mask_fg_tile)
1274 
1275  END SUBROUTINE adjust_nsst
1276 
1287  SUBROUTINE climo_trend(LATITUDE, MON, DAY, DELTSFC, DTREF)
1288  IMPLICIT NONE
1289 
1290  INTEGER, INTENT(IN) :: mon, day
1291 
1292  REAL, INTENT(IN) :: latitude, deltsfc
1293  REAL, INTENT(OUT) :: dtref
1294 
1295  INTEGER :: num_days(12), mon2, mon1
1296 
1297  REAL, TARGET :: sst_80_90(12)
1298  REAL, TARGET :: sst_70_80(12)
1299  REAL, TARGET :: sst_60_70(12)
1300  REAL, TARGET :: sst_50_60(12)
1301  REAL, TARGET :: sst_40_50(12)
1302  REAL, TARGET :: sst_30_40(12)
1303  REAL, TARGET :: sst_20_30(12)
1304  REAL, TARGET :: sst_10_20(12)
1305  REAL, TARGET :: sst_00_10(12)
1306  REAL, TARGET :: sst_m10_00(12)
1307  REAL, TARGET :: sst_m20_m10(12)
1308  REAL, TARGET :: sst_m30_m20(12)
1309  REAL, TARGET :: sst_m40_m30(12)
1310  REAL, TARGET :: sst_m50_m40(12)
1311  REAL, TARGET :: sst_m60_m50(12)
1312  REAL, TARGET :: sst_m70_m60(12)
1313  REAL, TARGET :: sst_m80_m70(12)
1314  REAL, TARGET :: sst_m90_m80(12)
1315 
1316  REAL, POINTER :: sst(:)
1317 
1318  DATA num_days /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
1319 
1320  DATA sst_80_90 /271.466, 271.458, 271.448, 271.445, 271.519, 271.636, &
1321  272.023, 272.066, 272.001, 271.698, 271.510, 271.472/
1322 
1323  DATA sst_70_80 /272.149, 272.103, 272.095, 272.126, 272.360, 272.988, &
1324  274.061, 274.868, 274.415, 273.201, 272.468, 272.268/
1325 
1326  DATA sst_60_70 /274.240, 274.019, 273.988, 274.185, 275.104, 276.875, &
1327  279.005, 280.172, 279.396, 277.586, 275.818, 274.803/
1328 
1329  DATA sst_50_60 /277.277, 276.935, 277.021, 277.531, 279.100, 281.357, &
1330  283.735, 285.171, 284.399, 282.328, 279.918, 278.199/
1331 
1332  DATA sst_40_50 /281.321, 280.721, 280.850, 281.820, 283.958, 286.588, &
1333  289.195, 290.873, 290.014, 287.652, 284.898, 282.735/
1334 
1335  DATA sst_30_40 /289.189, 288.519, 288.687, 289.648, 291.547, 293.904, &
1336  296.110, 297.319, 296.816, 295.225, 292.908, 290.743/
1337 
1338  DATA sst_20_30 /294.807, 294.348, 294.710, 295.714, 297.224, 298.703, &
1339  299.682, 300.127, 300.099, 299.455, 297.953, 296.177/
1340 
1341  DATA sst_10_20 /298.878, 298.720, 299.033, 299.707, 300.431, 300.709, &
1342  300.814, 300.976, 301.174, 301.145, 300.587, 299.694/
1343 
1344  DATA sst_00_10 /300.415, 300.548, 300.939, 301.365, 301.505, 301.141, &
1345  300.779, 300.660, 300.818, 300.994, 300.941, 300.675/
1346 
1347  DATA sst_m10_00 /300.226, 300.558, 300.914, 301.047, 300.645, 299.870, &
1348  299.114, 298.751, 298.875, 299.294, 299.721, 299.989/
1349 
1350  DATA sst_m20_m10 /299.547, 299.985, 300.056, 299.676, 298.841, 297.788, &
1351  296.893, 296.491, 296.687, 297.355, 298.220, 298.964/
1352 
1353  DATA sst_m30_m20 /297.524, 298.073, 297.897, 297.088, 295.846, 294.520, &
1354  293.525, 293.087, 293.217, 293.951, 295.047, 296.363/
1355 
1356  DATA sst_m40_m30 /293.054, 293.765, 293.468, 292.447, 291.128, 289.781, &
1357  288.773, 288.239, 288.203, 288.794, 289.947, 291.553/
1358 
1359  DATA sst_m50_m40 /285.052, 285.599, 285.426, 284.681, 283.761, 282.826, &
1360  282.138, 281.730, 281.659, 281.965, 282.768, 283.961/
1361 
1362  DATA sst_m60_m50 /277.818, 278.174, 277.991, 277.455, 276.824, 276.229, &
1363  275.817, 275.585, 275.560, 275.687, 276.142, 276.968/
1364 
1365  DATA sst_m70_m60 /273.436, 273.793, 273.451, 272.813, 272.349, 272.048, &
1366  271.901, 271.838, 271.845, 271.889, 272.080, 272.607/
1367 
1368  DATA sst_m80_m70 /271.579, 271.578, 271.471, 271.407, 271.392, 271.391, &
1369  271.390, 271.391, 271.394, 271.401, 271.422, 271.486/
1370 
1371  DATA sst_m90_m80 /271.350, 271.350, 271.350, 271.350, 271.350, 271.350, &
1372  271.350, 271.350, 271.350, 271.350, 271.350, 271.350/
1373 
1374  nullify(sst)
1375  IF (latitude > 80.0) THEN
1376  sst => sst_80_90
1377  ELSEIF (latitude > 70.0) THEN
1378  sst => sst_70_80
1379  ELSEIF (latitude > 60.0) THEN
1380  sst => sst_60_70
1381  ELSEIF (latitude > 50.0) THEN
1382  sst => sst_50_60
1383  ELSEIF (latitude > 40.0) THEN
1384  sst => sst_40_50
1385  ELSEIF (latitude > 30.0) THEN
1386  sst => sst_30_40
1387  ELSEIF (latitude > 20.0) THEN
1388  sst => sst_20_30
1389  ELSEIF (latitude > 10.0) THEN
1390  sst => sst_10_20
1391  ELSEIF (latitude > 0.0) THEN
1392  sst => sst_00_10
1393  ELSEIF (latitude > -10.0) THEN
1394  sst => sst_m10_00
1395  ELSEIF (latitude > -20.0) THEN
1396  sst => sst_m20_m10
1397  ELSEIF (latitude > -30.0) THEN
1398  sst => sst_m30_m20
1399  ELSEIF (latitude > -40.0) THEN
1400  sst => sst_m40_m30
1401  ELSEIF (latitude > -50.0) THEN
1402  sst => sst_m50_m40
1403  ELSEIF (latitude > -60.0) THEN
1404  sst => sst_m60_m50
1405  ELSEIF (latitude > -70.0) THEN
1406  sst => sst_m70_m60
1407  ELSEIF (latitude > -80.0) THEN
1408  sst => sst_m80_m70
1409  ELSE
1410  sst => sst_m90_m80
1411  END IF
1412 
1413  IF (day >= 15) THEN
1414  mon2 = mon+1
1415  IF(mon2 == 13) mon2 = 1
1416  mon1 = mon
1417  dtref = (sst(mon2) - sst(mon1)) / num_days(mon1)
1418  ELSE
1419  mon1 = mon - 1
1420  IF (mon1 == 0) mon1=12
1421  mon2 = mon
1422  dtref = (sst(mon2) - sst(mon1)) / num_days(mon1)
1423  ENDIF
1424 
1425  dtref = dtref * (deltsfc / 24.0)
1426 
1427  END SUBROUTINE climo_trend
1428 
1440  SUBROUTINE dtzm_point(XT,XZ,DT_COOL,ZC,Z1,Z2,DTZM)
1441  implicit none
1442 
1443  real, intent(in) :: xt,xz,dt_cool,zc,z1,z2
1444  real, intent(out) :: dtzm
1445 
1446  real, parameter :: zero = 0.0
1447  real, parameter :: one = 1.0
1448  real, parameter :: half = 0.5
1449  real :: dt_warm,dtw,dtc
1450 !
1451 ! get the mean warming in the range of z=z1 to z=z2
1452 !
1453  dtw = zero
1454  if ( xt > zero ) then
1455  dt_warm = (xt+xt)/xz ! Tw(0)
1456  if ( z1 < z2) then
1457  if ( z2 < xz ) then
1458  dtw = dt_warm*(one-(z1+z2)/(xz+xz))
1459  elseif ( z1 < xz .and. z2 >= xz ) then
1460  dtw = half*(one-z1/xz)*dt_warm*(xz-z1)/(z2-z1)
1461  endif
1462  elseif ( z1 == z2 ) then
1463  if ( z1 < xz ) then
1464  dtw = dt_warm*(one-z1/xz)
1465  endif
1466  endif
1467  endif
1468 !
1469 ! get the mean cooling in the range of z=z1 to z=z2
1470 !
1471  dtc = zero
1472  if ( zc > zero ) then
1473  if ( z1 < z2) then
1474  if ( z2 < zc ) then
1475  dtc = dt_cool*(one-(z1+z2)/(zc+zc))
1476  elseif ( z1 < zc .and. z2 >= zc ) then
1477  dtc = half*(one-z1/zc)*dt_cool*(zc-z1)/(z2-z1)
1478  endif
1479  elseif ( z1 == z2 ) then
1480  if ( z1 < zc ) then
1481  dtc = dt_cool*(one-z1/zc)
1482  endif
1483  endif
1484  endif
1485 
1486 !
1487 ! get the mean T departure from Tf in the range of z=z1 to z=z2
1488 !
1489  dtzm = dtw - dtc
1490 
1491  END SUBROUTINE dtzm_point
1492 
1512  subroutine tf_thaw_set(tf_ij,mask_ij,itile,jtile,tice,tclm,tf_thaw,nx,ny, &
1513  nset_thaw_s,nset_thaw_i,nset_thaw_c)
1514 
1515  real, dimension(nx*ny), intent(in) :: tf_ij
1516  integer, dimension(nx*ny), intent(in) :: mask_ij
1517  real, intent(in) :: tice,tclm
1518  integer, intent(in) :: itile,jtile,nx,ny
1519  real, intent(out) :: tf_thaw
1520  integer, intent(inout) :: nset_thaw_s,nset_thaw_i,nset_thaw_c
1521 ! Local
1522  real, parameter :: bmiss = -999.0
1523  real, dimension(nx,ny) :: tf
1524  integer, dimension(nx,ny) :: mask
1525  integer :: krad,max_search,istart,iend,jstart,jend
1526  integer :: ii,jj,iii,jjj
1527  logical :: is_ice
1528 
1529  max_search = 2
1530 
1531  mask(:,:) = reshape(mask_ij,(/nx,ny/) )
1532  tf(:,:) = reshape(tf_ij,(/nx,ny/) )
1533 
1534  tf_thaw = bmiss
1535 
1536  do krad = 1, max_search
1537 
1538  istart = itile - krad
1539  iend = itile + krad
1540  jstart = jtile - krad
1541  jend = jtile + krad
1542 
1543  do jj = jstart, jend
1544  do ii = istart, iend
1545 
1546 
1547  if ((jj == jstart) .or. (jj == jend) .or. &
1548  (ii == istart) .or. (ii == iend)) then
1549 
1550  if ((jj >= 1) .and. (jj <= ny)) then
1551  jjj = jj
1552  if (ii <= 0) then
1553  iii = nx + ii
1554  else if (ii >= (nx+1)) then
1555  iii = ii - nx
1556  else
1557  iii = ii
1558  endif
1559 
1560 !----------------------------------------------------------------------
1561 ! SEE IF NEARBY POINTS ARE SEA ICE. IF THEY ARE, AND THE SEARCH FOR
1562 ! A GAUSSIAN GRID OPEN WATER POINT FAILS, THEN TREF WILL BE SET TO
1563 ! FREEZING BELOW.
1564 !----------------------------------------------------------------------
1565  if (krad <= 2 .and. mask(iii,jjj) == 2) is_ice = .true.
1566 
1567  if (mask(iii,jjj) == 0) then
1568  tf_thaw = tf(iii,jjj)
1569  nset_thaw_s = nset_thaw_s + 1
1570  write(*,'(a,I4,2F9.3)') 'nset_thaw_s,tf(iii,jjj),tclm : ',nset_thaw_s,tf(iii,jjj),tclm
1571  go to 100
1572  endif ! tile mask is open water
1573 
1574  endif
1575  endif
1576  enddo
1577  enddo
1578  enddo ! krad loop
1579 
1580  100 continue
1581 
1582  if ( tf_thaw == bmiss ) then
1583  if (is_ice) then
1584  tf_thaw = tice
1585  nset_thaw_i = nset_thaw_i + 1
1586  write(*,'(a,I4,F9.3)') 'nset_thaw_i,tf_ice : ',nset_thaw_i,tice
1587  else
1588  tf_thaw = 0.8*tice+0.2*tclm
1589  nset_thaw_c = nset_thaw_c + 1
1590  write(*,'(a,I4,2F9.3)') 'nset_thaw_c,tf_ice,tclm : ',nset_thaw_c,tice,tclm
1591  endif
1592  endif
1593 
1594  end subroutine tf_thaw_set
1595 
1603  subroutine nsst_water_reset(nsst,ij,tf_thaw)
1604  use read_write_data, only : nsst_data
1605  implicit none
1606 
1607  integer, intent(in) :: ij
1608 
1609  real, intent(in) :: tf_thaw
1610 
1611  type(nsst_data), intent(inout) :: nsst
1612 
1613  nsst%c_0(ij) = 0.0
1614  nsst%c_d(ij) = 0.0
1615  nsst%d_conv(ij) = 0.0
1616  nsst%dt_cool(ij) = 0.0
1617  nsst%ifd(ij) = 0.0
1618  nsst%qrain(ij) = 0.0
1619  nsst%tref(ij) = tf_thaw
1620  nsst%w_0(ij) = 0.0
1621  nsst%w_d(ij) = 0.0
1622  nsst%xs(ij) = 0.0
1623  nsst%xt(ij) = 0.0
1624  nsst%xtts(ij) = 0.0
1625  nsst%xu(ij) = 0.0
1626  nsst%xv(ij) = 0.0
1627  nsst%xz(ij) = 30.0
1628  nsst%xzts(ij) = 0.0
1629  nsst%z_c(ij) = 0.0
1630  nsst%zm(ij) = 0.0
1631 
1632  end subroutine nsst_water_reset
1633 
1649 subroutine get_tf_clm(xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,tf_clm,tf_trd)
1650  use read_write_data, only : get_tf_clm_dim
1651 
1652  implicit none
1653 
1654  real, dimension(nx*ny), intent(in) :: xlats_ij
1655  real, dimension(nx*ny), intent(in) :: xlons_ij
1656  real, dimension(nx,ny), intent(out) :: tf_clm
1657  real, dimension(nx,ny), intent(out) :: tf_trd
1658  integer, intent(in) :: iy,im,id,ih,nx,ny
1659 ! local declare
1660  real, allocatable, dimension(:,:) :: tf_clm0 ! sst climatology at the valid time (nxc,nyc)
1661  real, allocatable, dimension(:,:) :: tf_trd0 ! 6-hourly sst climatology tendency valid at atime (nxc,nyc)
1662  real, allocatable, dimension(:) :: cxlats ! latitudes of sst climatology
1663  real, allocatable, dimension(:) :: cxlons ! longitudes of sst climatology
1664 
1665  real, dimension(nx*ny) :: tf_clm_ij ! sst climatology at target grids (nx*ny)
1666  real, dimension(nx*ny) :: tf_trd_ij ! 6-hourly sst climatology tendency
1667  real :: wei1,wei2
1668  integer :: nxc,nyc,mon1,mon2
1669  character (len=6), parameter :: fin_tf_clm='sstclm' ! sst climatology file name
1670 !
1671 ! get which two months used and their weights from atime
1672 !
1673  call get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1674 !
1675 ! get the dimensions of the sst climatology & allocate the related arrays
1676 !
1677  call get_tf_clm_dim(fin_tf_clm,nyc,nxc)
1678  allocate( tf_clm0(nxc,nyc),tf_trd0(nxc,nyc),cxlats(nyc),cxlons(nxc) )
1679 !
1680 ! get tf_clm at the analysis time from monthly climatology & cxlats, cxlons
1681 !
1682  call get_tf_clm_ta(tf_clm0,tf_trd0,cxlats,cxlons,nyc,nxc,mon1,mon2,wei1,wei2)
1683 !
1684 ! get tf_clm (nx by ny lat/lon) valid at atime
1685 !
1686  if ( nx == nxc .and. ny == nyc ) then
1687  tf_clm(:,:) = tf_clm0(:,:)
1688  tf_trd(:,:) = tf_trd0(:,:)
1689 ! write(*,'(a,2f9.3)') 'same dimensions, tf_clm, min : ',minval(tf_clm),maxval(tf_clm)
1690  else
1691 ! write(*,'(a,4i8)') 'different dimensions,nx,ny,nxc,nyc : ',nx,ny,nxc,nyc
1692  call intp_tile(tf_clm0, cxlats, cxlons, nyc, nxc, &
1693  tf_clm_ij,xlats_ij,xlons_ij,ny, nx)
1694  call intp_tile(tf_trd0, cxlats, cxlons, nyc, nxc, &
1695  tf_trd_ij,xlats_ij,xlons_ij,ny, nx)
1696 ! write(*,'(a,2f9.3)') 'tf_clm0, min, max : ',minval(tf_clm0),maxval(tf_clm0)
1697 
1698  tf_clm(:,:) = reshape(tf_clm_ij, (/nx,ny/) )
1699  tf_trd(:,:) = reshape(tf_trd_ij, (/nx,ny/) )
1700  endif
1701 
1702 end subroutine get_tf_clm
1703 
1718 subroutine get_tf_clm_ta(tf_clm_ta,tf_clm_trend,xlats,xlons,nlat,nlon,mon1,mon2,wei1,wei2)
1719  use read_write_data, only : read_tf_clim_grb
1720  implicit none
1721 
1722 ! input
1723  integer, intent(in) :: nlat,nlon,mon1,mon2
1724  real, intent(in) :: wei1,wei2
1725 ! output
1726  real, dimension(nlon,nlat), intent(out) :: tf_clm_ta,tf_clm_trend
1727  real, dimension(nlat), intent(out) :: xlats
1728  real, dimension(nlon), intent(out) :: xlons
1729 
1730 !input/output data file names
1731  character (len=6), parameter :: fin_tf_clm='sstclm'
1732 
1733 ! local declare
1734  real, dimension(nlon,nlat) :: tf_clm1,tf_clm2
1735 
1736 !
1737 ! read in rtg sst climatology without bitmap (surface mask) for mon1 and mon2
1738 !
1739  call read_tf_clim_grb(trim(fin_tf_clm),tf_clm1,xlats,xlons,nlat,nlon,mon1)
1740  call read_tf_clim_grb(trim(fin_tf_clm),tf_clm2,xlats,xlons,nlat,nlon,mon2)
1741 !
1742 ! tf_clim at the analysis time
1743 !
1744  tf_clm_ta(:,:) = wei1*tf_clm1(:,:)+wei2*tf_clm2(:,:)
1745 !
1746 ! tf_clim trend at the analysis time (month to month)
1747 !
1748  tf_clm_trend(:,:) = (tf_clm2(:,:)-tf_clm1(:,:))/120.0
1749 
1750  write(*,'(a,2f9.3)') 'tf_clm_ta, min, max : ',minval(tf_clm_ta),maxval(tf_clm_ta)
1751  write(*,'(a,2f9.3)') 'tf_clm_trend, min, max : ',minval(tf_clm_trend),maxval(tf_clm_trend)
1752  end subroutine get_tf_clm_ta
1753 
1766 subroutine get_sal_clm(xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,sal_clm)
1767  use read_write_data, only : get_dim_nc
1768  implicit none
1769 
1770  real, dimension(nx*ny), intent(in) :: xlats_ij !
1771  real, dimension(nx*ny), intent(in) :: xlons_ij !
1772  real, dimension(nx,ny), intent(out) :: sal_clm !
1773  integer, intent(in) :: iy,im,id,ih,nx,ny
1774 ! local declare
1775  real, allocatable, dimension(:,:) :: sal_clm0 ! salinity climatology at the valid time
1776  real, allocatable, dimension(:) :: cxlats ! latitudes of sst climatology
1777  real, allocatable, dimension(:) :: cxlons ! longitudes of sst climatology
1778 
1779  real, dimension(nx*ny) :: sal_clm_ij ! salinity climatology at target grids (nx*ny)
1780  real :: wei1,wei2
1781  integer :: nxc,nyc,mon1,mon2
1782  character (len=6), parameter :: fin_sal_clm='salclm' ! salinity climatology file name
1783 !
1784 ! get which two months used and their weights from atime
1785 !
1786  call get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1787 !
1788 ! get the dimensions of the sst climatology & allocate the related arrays
1789 !
1790  call get_dim_nc(fin_sal_clm,nyc,nxc)
1791  allocate( sal_clm0(nxc,nyc),cxlats(nyc),cxlons(nxc) )
1792 !
1793 ! get sal_clm at the analysis time from monthly climatology & cxlats, cxlons
1794 !
1795  call get_sal_clm_ta(sal_clm0,cxlats,cxlons,nyc,nxc,mon1,mon2,wei1,wei2)
1796 !
1797 ! get sal_clm (nx by ny lat/lon) valid at atime
1798 !
1799  if ( nx == nxc .and. ny == nyc ) then
1800  sal_clm(:,:) = sal_clm0(:,:)
1801 ! write(*,'(a,2f9.3)') 'same dimensions, sal_clm, min : ',minval(sal_clm),maxval(sal_clm)
1802  else
1803 ! write(*,'(a,4i8)') 'different dimensions,nx,ny,nxc,nyc : ',nx,ny,nxc,nyc
1804  call intp_tile(sal_clm0, cxlats, cxlons, nyc, nxc, &
1805  sal_clm_ij,xlats_ij,xlons_ij,ny, nx)
1806 ! write(*,'(a,2f9.3)') 'sal_clm0, min, max : ',minval(sal_clm0),maxval(sal_clm0)
1807 ! write(*,'(a,2f9.3)') 'done with intp_tile for sal_clm, min, max : ',minval(sal_clm_ij),maxval(sal_clm_ij)
1808 
1809  sal_clm(:,:) = reshape(sal_clm_ij, (/nx,ny/) )
1810  endif
1811 
1812 end subroutine get_sal_clm
1813 
1826 subroutine get_sal_clm_ta(sal_clm_ta,xlats,xlons,nlat,nlon,mon1,mon2,wei1,wei2)
1827 
1828  use read_write_data, only : read_salclm_gfs_nc
1829  implicit none
1830 
1831 ! input
1832  integer, intent(in) :: nlat,nlon,mon1,mon2
1833  real, intent(in) :: wei1,wei2
1834 ! output
1835  real, dimension(nlon,nlat), intent(out) :: sal_clm_ta
1836  real, dimension(nlat), intent(out) :: xlats
1837  real, dimension(nlon), intent(out) :: xlons
1838 
1839 !input/output data file names
1840  character (len=6), parameter :: fin_sal_clm='salclm'
1841 
1842 ! local declare
1843  real, dimension(nlon,nlat) :: sal_clm1,sal_clm2
1844 
1845 !
1846 ! read in rtg sst climatology without bitmap (surface mask) for mon1 and mon2
1847 !
1848  call read_salclm_gfs_nc(trim(fin_sal_clm),sal_clm1,xlats,xlons,nlat,nlon,mon1)
1849  call read_salclm_gfs_nc(trim(fin_sal_clm),sal_clm2,xlats,xlons,nlat,nlon,mon2)
1850 !
1851 ! sal_clim at the analysis time
1852 !
1853  sal_clm_ta(:,:) = wei1*sal_clm1(:,:)+wei2*sal_clm2(:,:)
1854  write(*,'(a,2f9.3)') 'sal_clm_ta, min, max : ',minval(sal_clm_ta),maxval(sal_clm_ta)
1855  end subroutine get_sal_clm_ta
1856 
1871 subroutine intp_tile(tf_lalo,dlats_lalo,dlons_lalo,jdim_lalo,idim_lalo, &
1872  tf_tile,xlats_tile,xlons_tile,jdim_tile,idim_tile)
1873 
1874  use utils
1875 
1876  implicit none
1877 
1878 ! input/output
1879  real, dimension(idim_lalo,jdim_lalo), intent(in) :: tf_lalo
1880  real, dimension(jdim_lalo), intent(in) :: dlats_lalo
1881  real, dimension(idim_lalo), intent(in) :: dlons_lalo
1882  real, dimension(jdim_tile*idim_tile), intent(in) :: xlats_tile
1883  real, dimension(jdim_tile*idim_tile), intent(in) :: xlons_tile
1884  integer, intent(in) :: jdim_lalo,idim_lalo,jdim_tile,idim_tile
1885  real, dimension(jdim_tile*idim_tile), intent(out) :: tf_tile
1886 
1887 ! local
1888  real, parameter :: deg2rad=3.1415926/180.0
1889  real, dimension(jdim_lalo) :: xlats_lalo
1890  real, dimension(idim_lalo) :: xlons_lalo
1891  real :: wsum
1892  integer :: itile,jtile
1893  integer :: ij
1894  integer :: ilalo,jlalo,ilalop1,jlalop1
1895 
1896  integer, allocatable, dimension(:,:) :: id1,id2,jdc
1897  real, allocatable, dimension(:,:,:) :: agrid,s2c
1898 
1899  print*
1900  print*,'interpolate from lat/lon grids to any one grid with known lat/lon'
1901 
1902  xlats_lalo = dlats_lalo*deg2rad
1903  xlons_lalo = dlons_lalo*deg2rad
1904 
1905  allocate(agrid(idim_tile,jdim_tile,2))
1906  agrid(:,:,1) = reshape(xlons_tile, (/idim_tile,jdim_tile/) )
1907  agrid(:,:,2) = reshape(xlats_tile, (/idim_tile,jdim_tile/) )
1908  agrid = agrid*deg2rad
1909 
1910  allocate(id1(idim_tile,jdim_tile))
1911  allocate(id2(idim_tile,jdim_tile))
1912  allocate(jdc(idim_tile,jdim_tile))
1913  allocate(s2c(idim_tile,jdim_tile,4))
1914 
1915 !----------------------------------------------------------------------
1916 ! compute bilinear weights for each model point from the nearest
1917 ! four lalo points. does not account for mask. that
1918 ! happens later.
1919 !----------------------------------------------------------------------
1920 
1921  call remap_coef( 1, idim_tile, 1, jdim_tile, idim_lalo, jdim_lalo, &
1922  xlons_lalo, xlats_lalo, id1, id2, jdc, s2c, agrid )
1923 
1924  do ij = 1, jdim_tile*idim_tile
1925 
1926  jtile = (ij-1)/idim_tile + 1
1927  itile = mod(ij,idim_tile)
1928  if (itile==0) itile = idim_tile
1929 
1930  ilalo = id1(itile,jtile)
1931  ilalop1 = id2(itile,jtile)
1932  jlalo = jdc(itile,jtile)
1933  jlalop1 = jdc(itile,jtile) + 1
1934 
1935  wsum = s2c(itile,jtile,1) + s2c(itile,jtile,2) + &
1936  s2c(itile,jtile,3) + s2c(itile,jtile,4)
1937 
1938  tf_tile(ij) = ( s2c(itile,jtile,1)*tf_lalo(ilalo,jlalo) + &
1939  s2c(itile,jtile,2)*tf_lalo(ilalop1,jlalo) + &
1940  s2c(itile,jtile,3)*tf_lalo(ilalop1,jlalop1) + &
1941  s2c(itile,jtile,4)*tf_lalo(ilalo,jlalop1) )/wsum
1942  enddo
1943 
1944  deallocate(id1, id2, jdc, s2c)
1945 
1946 end subroutine intp_tile
1947 
1960 subroutine get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1961  implicit none
1962 
1963 ! input
1964  integer, intent(in) :: iy,im,id,ih
1965 ! output
1966  integer, intent(out) :: mon1,mon2
1967  real, intent(out) :: wei1,wei2
1968 
1969 ! local declare
1970  real :: rjday
1971  integer :: mon,monend,monm,monp,jdow,jdoy,jday
1972  integer :: jda(8)
1973 !
1974 !dayhf : julian day of the middle of each month
1975 !
1976  real, dimension(13) :: dayhf
1977  data dayhf/15.5,45.0,74.5,105.0,135.5,166.0,196.5,227.5,258.0,288.5,319.0,349.5,380.5/
1978 
1979 ! 15, 44, 73.5, 104, 134.5, 165, 195.5, 226.5, 257, 287.5, 318.5, 349 ! from
1980 ! woa05
1981 
1982  jda=0
1983  jda(1)=iy
1984  jda(2)=im
1985  jda(3)=id
1986  jda(5)=ih
1987  jdow = 0
1988  jdoy = 0
1989  jday = 0
1990  call w3doxdat(jda,jdow,jdoy,jday)
1991  rjday=jdoy+jda(5)/24.
1992  if(rjday.lt.dayhf(1)) rjday=rjday+365.
1993 
1994  monend = 12
1995  do mon = 1, monend
1996  monm = mon
1997  monp = mon + 1
1998  if( rjday >= dayhf(monm) .and. rjday < dayhf(monp) ) then
1999  mon1 = monm
2000  mon2 = monp
2001  go to 10
2002  endif
2003  enddo
2004 
2005  print *,'FATAL ERROR in get_tim_wei, wrong rjday',rjday
2006  call abort
2007  10 continue
2008 
2009  wei1 = (dayhf(mon2)-rjday)/(dayhf(mon2)-dayhf(mon1))
2010  wei2 = (rjday-dayhf(mon1))/(dayhf(mon2)-dayhf(mon1))
2011 
2012  if( mon2 == 13 ) mon2=1
2013 
2014  write(*,'(a,2i4,3f9.3)') 'mon1,mon2,rjday,wei1,wei2=',mon1,mon2,rjday,wei1,wei2
2015 
2016  end subroutine get_tim_wei
2017 
2027  real function tfreez(salinity)
2028 
2029  implicit none
2030 
2031  REAL salinity,sal
2032  REAL a1, a2, a3
2033  parameter(a1 = -0.0575)
2034  parameter(a2 = 1.710523e-3)
2035  parameter(a3 = -2.154996e-4)
2036 
2037  IF (salinity .LT. 0.) THEN
2038  sal = 0.
2039  ELSE
2040  sal = salinity
2041  ENDIF
2042  tfreez = sal*(a1+a2*sqrt(sal)+a3*sal)
2043 
2044  return
2045  end
num_parthds
integer function num_parthds()
Return the number of omp threads.
Definition: num_parthds.f90:9
land_increments::add_increment_soil
subroutine, public add_increment_soil(rla, rlo, stc_state, smc_state, slc_state, stc_updated, slc_updated, soilsnow_tile, soilsnow_fg_tile, lensfc, lsoil, idim, jdim, lsm, myrank)
Read in gsi file with soil state increments (on the gaussian grid), interpolate increments to the cub...
Definition: land_increments.f90:55
get_tf_clm
subroutine get_tf_clm(xlats_ij, xlons_ij, ny, nx, iy, im, id, ih, tf_clm, tf_trd)
Get the sst climatology at the valid time and on the target grid.
Definition: cycle.f90:1649
read_write_data::read_lat_lon_orog
subroutine, public read_lat_lon_orog(RLA, RLO, OROG, OROG_UF, TILE_NUM, IDIM, JDIM, IJDIM, LANDFRAC)
Read latitude and longitude for the cubed-sphere tile from the &#39;grid&#39; file.
Definition: read_write_data.f90:680
get_tf_clm_ta
subroutine get_tf_clm_ta(tf_clm_ta, tf_clm_trend, xlats, xlons, nlat, nlon, mon1, mon2, wei1, wei2)
Get the reference temperature/sst climatology and its trend at analysis time.
Definition: cycle.f90:1718
dtzm_point
subroutine dtzm_point(XT, XZ, DT_COOL, ZC, Z1, Z2, DTZM)
Compute the vertical mean of the NSST t-profile.
Definition: cycle.f90:1440
get_sal_clm_ta
subroutine get_sal_clm_ta(sal_clm_ta, xlats, xlons, nlat, nlon, mon1, mon2, wei1, wei2)
Get climatological salinity at the analysis time.
Definition: cycle.f90:1826
sfcdrv
subroutine sfcdrv(LUGB, IDIM, JDIM, LENSFC, LSOIL, DELTSFC, IY, IM, ID, IH, FH, IALB, USE_UFO, DO_NSST, DO_SFCCYCLE, DO_LNDINC, FRAC_GRID, ZSEA1, ZSEA2, ISOT, IVEGSRC, MYRANK)
Driver routine for updating the surface fields.
Definition: cycle.f90:299
tf_thaw_set
subroutine tf_thaw_set(tf_ij, mask_ij, itile, jtile, tice, tclm, tf_thaw, nx, ny, nset_thaw_s, nset_thaw_i, nset_thaw_c)
Set the background reference temperature (tf) for points where the ice has just melted.
Definition: cycle.f90:1512
land_increments::add_increment_snow
subroutine, public add_increment_snow(snd_inc, mask, lensfc, snd)
Add snow depth increment to model snow depth state, and limit output to be non-negative.
Definition: land_increments.f90:378
land_increments::calculate_landinc_mask
subroutine, public calculate_landinc_mask(smc, swe, vtype, lensfc, veg_type_landice, mask)
Calculate soil mask for land on model grid.
Definition: land_increments.f90:408
land_increments::apply_land_da_adjustments_snd
subroutine, public apply_land_da_adjustments_snd(lsm, lensfc, mask, swe_bck, snd_bck, snd_anl, swe_adj)
Make adjustments to dependent variables after applying land increments.
Definition: land_increments.f90:626
read_write_data::get_dim_nc
subroutine, public get_dim_nc(filename, nlat, nlon)
Get the i/j dimensions of the data from a NetCDF file.
Definition: read_write_data.f90:1824
machine
Holds machine dependent constants for global_cycle.
Definition: machine.f90:7
read_write_data::read_salclm_gfs_nc
subroutine, public read_salclm_gfs_nc(filename, sal, xlats, xlons, nlat, nlon, itime)
Read the woa05 salinity monthly climatology file.
Definition: read_write_data.f90:1753
land_increments
Definition: land_increments.f90:5
read_write_data::nsst_data
Definition: read_write_data.f90:16
intp_tile
subroutine intp_tile(tf_lalo, dlats_lalo, dlons_lalo, jdim_lalo, idim_lalo, tf_tile, xlats_tile, xlons_tile, jdim_tile, idim_tile)
Interpolate lon/lat grid data to the fv3 native grid (tf_lalo =&gt; tf_tile).
Definition: cycle.f90:1871
land_increments::apply_land_da_adjustments_soil
subroutine, public apply_land_da_adjustments_soil(lsm, isot, ivegsrc, lensfc, lsoil, rsoiltype, mask, stc_bck, stc_adj, smc_adj, slc_adj, stc_updated, slc_updated, zsoil)
Make adjustments to dependent variables after applying land increments.
Definition: land_increments.f90:465
climo_trend
subroutine climo_trend(LATITUDE, MON, DAY, DELTSFC, DTREF)
If the tile point is an isolated water point that has no corresponding gsi water point, then tref is updated using the rtg sst climo trend.
Definition: cycle.f90:1287
read_write_data::read_data
subroutine, public read_data(LSOIL, LENSFC, DO_NSST, INC_FILE, IS_NOAHMP, TSFFCS, SMCFCS, SWEFCS, STCFCS, TG3FCS, ZORFCS, CVFCS, CVBFCS, CVTFCS, ALBFCS, VEGFCS, SLIFCS, CNPFCS, F10M, VETFCS, SOTFCS, ALFFCS, USTAR, FMM, FHH, SIHFCS, SICFCS, SITFCS, TPRCP, SRFLAG, SNDFCS, VMNFCS, VMXFCS, SLCFCS, SLPFCS, ABSFCS, T2M, Q2M, SLMASK, ZSOIL, NSST)
Read the first guess surface records and nsst records (if selected) for a single cubed-sphere tile...
Definition: read_write_data.f90:1029
utils::remap_coef
subroutine, public remap_coef(is, ie, js, je, im, jm, lon, lat, id1, id2, jdc, s2c, agrid)
Generate the weights and index of the grids used in the bilinear interpolation.
Definition: utils.F90:38
read_write_data::get_tf_clm_dim
subroutine, public get_tf_clm_dim(file_sst, mlat_sst, mlon_sst)
Get the i/j dimensions of RTG SST climatology file.
Definition: read_write_data.f90:1695
sfc_drv
program sfc_drv
Stand alone surface/NSST cycle driver for the cubed-sphere grid.
Definition: cycle.f90:101
adjust_nsst
subroutine adjust_nsst(RLA, RLO, SLMSK_TILE, SLMSK_FG_TILE, SKINT_TILE, SICET_TILE, sice_tile, sice_fg_tile, SOILT_TILE, NSST, LENSFC, LSOIL, IDIM, JDIM, ZSEA1, ZSEA2, tf_clm_tile, tf_trd_tile, sal_clm_tile, LANDFRAC, FRAC_GRID)
Read in gsi file with the updated reference temperature increments (on the gaussian grid)...
Definition: cycle.f90:859
read_write_data::read_gsi_data
subroutine, public read_gsi_data(GSI_FILE, FILE_TYPE, LSOIL)
Read increment file from the GSI containing either the foundation temperature increments and mask...
Definition: read_write_data.f90:857
tfreez
real function tfreez(salinity)
Compute the freezing point of water as a function of salinity.
Definition: cycle.f90:2027
read_write_data
This module contains routines that read and write data.
Definition: read_write_data.f90:8
get_tim_wei
subroutine get_tim_wei(iy, im, id, ih, mon1, mon2, wei1, wei2)
For a given date, determine the bounding months and the linear time interpolation weights...
Definition: cycle.f90:1960
get_sal_clm
subroutine get_sal_clm(xlats_ij, xlons_ij, ny, nx, iy, im, id, ih, sal_clm)
Get salinity climatology at the valid time on the target grid.
Definition: cycle.f90:1766
nsst_water_reset
subroutine nsst_water_reset(nsst, ij, tf_thaw)
If the first guess was sea ice, but the analysis is open water, reset all nsst variables.
Definition: cycle.f90:1603
read_write_data::write_data
subroutine, public write_data(lensfc, idim, jdim, lsoil, do_nsst, nsst, slifcs, tsffcs, vegfcs, swefcs, tg3fcs, zorfcs, albfcs, alffcs, cnpfcs, f10m, t2m, q2m, vetfcs, sotfcs, ustar, fmm, fhh, sicfcs, sihfcs, sitfcs, tprcp, srflag, swdfcs, vmnfcs, vmxfcs, slpfcs, absfcs, slcfcs, smcfcs, stcfcs)
Update surface records - and nsst records if selected - on a single cubed-sphere tile to a pre-existi...
Definition: read_write_data.f90:120
utils
Module containing utility routines.
Definition: utils.F90:7
read_write_data::read_tf_clim_grb
subroutine, public read_tf_clim_grb(file_sst, sst, rlats_sst, rlons_sst, mlat_sst, mlon_sst, mon)
Read a GRIB1 sst climatological analysis file.
Definition: read_write_data.f90:1553