global_cycle  1.5.0
 All Data Structures Files Functions Variables
cycle.f90
Go to the documentation of this file.
1 
4 
89  PROGRAM sfc_drv
90 
91  use mpi
92 
93  IMPLICIT NONE
94 !
95  CHARACTER(LEN=3) :: donst
96  INTEGER :: idim, jdim, lsoil, lugb, iy, im, id, ih, ialb
97  INTEGER :: isot, ivegsrc, lensfc, zsea1_mm, zsea2_mm, ierr
98  INTEGER :: nprocs, myrank, num_threads, num_parthds, max_tasks
99  REAL :: fh, deltsfc, zsea1, zsea2
100  LOGICAL :: use_ufo, do_nsst, adjt_nst_only
101 !
102  namelist/namcyc/ idim,jdim,lsoil,lugb,iy,im,id,ih,fh, &
103  deltsfc,ialb,use_ufo,donst, &
104  adjt_nst_only,isot,ivegsrc,zsea1_mm, &
105  zsea2_mm, max_tasks
106 !
107  DATA idim,jdim,lsoil/96,96,4/
108  DATA iy,im,id,ih,fh/1997,8,2,0,0./
109  DATA lugb/51/, deltsfc/0.0/, ialb/1/, max_tasks/99999/
110  DATA isot/1/, ivegsrc/2/, zsea1_mm/0/, zsea2_mm/0/
111 !
112  CALL mpi_init(ierr)
113  CALL mpi_comm_size(mpi_comm_world, nprocs, ierr)
114  CALL mpi_comm_rank(mpi_comm_world, myrank, ierr)
115 
116  if (myrank==0) call w3tagb('GLOBAL_CYCLE',2018,0179,0055,'NP20')
117 
118  num_threads = num_parthds()
119 
120  print*
121  print*,"STARTING CYCLE PROGRAM ON RANK ", myrank
122  print*,"RUNNING WITH ", nprocs, "TASKS"
123  print*,"AND WITH ", num_threads, " THREADS."
124 
125  use_ufo = .false.
126  donst = "NO"
127  adjt_nst_only = .false.
128 
129  print*
130  print*,"READ NAMCYC NAMELIST."
131 
132  CALL baopenr(36, "fort.36", ierr)
133  READ(36, nml=namcyc)
134  IF (myrank==0) WRITE(6,namcyc)
135 
136  IF (max_tasks < 99999 .AND. myrank > (max_tasks - 1)) THEN
137  print*,"USER SPECIFIED MAX NUMBER OF TASKS: ", max_tasks
138  print*,"WILL NOT RUN CYCLE PROGRAM ON RANK: ", myrank
139  goto 333
140  ENDIF
141 
142  lensfc = idim*jdim ! TOTAL NUMBER OF POINTS FOR THE CUBED-SPHERE TILE
143 
144  zsea1 = float(zsea1_mm) / 1000.0 ! CONVERT FROM MM TO METERS
145  zsea2 = float(zsea2_mm) / 1000.0
146 
147  IF (donst == "YES") THEN
148  do_nsst=.true.
149  ELSE
150  do_nsst=.false.
151  ENDIF
152 
153  print*
154  IF (myrank==0) print*,"LUGB,IDIM,JDIM,LSOIL,DELTSFC,IY,IM,ID,IH,FH: ", &
155  lugb,idim,jdim,lsoil,deltsfc,iy,im,id,ih,fh
156 
157  CALL sfcdrv(lugb,idim,jdim,lensfc,lsoil,deltsfc, &
158  iy,im,id,ih,fh,ialb, &
159  use_ufo,do_nsst,adjt_nst_only, &
160  zsea1,zsea2,isot,ivegsrc,myrank)
161 
162  print*
163  print*,'CYCLE PROGRAM COMPLETED NORMALLY ON RANK: ', myrank
164 
165  333 CONTINUE
166 
167  CALL mpi_barrier(mpi_comm_world, ierr)
168 
169  if (myrank==0) call w3tage('GLOBAL_CYCLE')
170 
171  CALL mpi_finalize(ierr)
172 
173  stop
174 
175  END PROGRAM sfc_drv
176 
283  SUBROUTINE sfcdrv(LUGB, IDIM,JDIM,LENSFC,LSOIL,DELTSFC, &
284  iy,im,id,ih,fh,ialb, &
285  use_ufo,do_nsst,adjt_nst_only, &
286  zsea1,zsea2,isot,ivegsrc,myrank)
287 !
288  USE read_write_data
289 
290  IMPLICIT NONE
291 
292  INTEGER, INTENT(IN) :: idim, jdim, lensfc, lsoil, ialb
293  INTEGER, INTENT(IN) :: lugb, iy, im, id, ih
294  INTEGER, INTENT(IN) :: isot, ivegsrc, myrank
295 
296  LOGICAL, INTENT(IN) :: use_ufo, do_nsst, adjt_nst_only
297 
298  REAL, INTENT(IN) :: fh, deltsfc, zsea1, zsea2
299 
300  INTEGER, PARAMETER :: nlunit=35
301  INTEGER, PARAMETER :: sz_nml=1
302 
303  CHARACTER(LEN=5) :: tile_num
304  CHARACTER(LEN=500) :: gsi_file
305  CHARACTER(LEN=4) :: input_nml_file(sz_nml)
306 
307  INTEGER :: i, ierr
308  INTEGER :: i_index(lensfc), j_index(lensfc)
309  INTEGER :: idum(idim,jdim)
310 
311  REAL :: slmask(lensfc), orog(lensfc)
312  REAL :: sihfcs(lensfc), sicfcs(lensfc)
313  REAL :: sitfcs(lensfc), tsffcs(lensfc)
314  REAL :: snofcs(lensfc), zorfcs(lensfc)
315  REAL :: albfcs(lensfc,4), tg3fcs(lensfc)
316  REAL :: cnpfcs(lensfc), smcfcs(lensfc,lsoil)
317  REAL :: stcfcs(lensfc,lsoil), slifcs(lensfc)
318  REAL :: aisfcs(lensfc), f10m(lensfc)
319  REAL :: vegfcs(lensfc), vetfcs(lensfc)
320  REAL :: sotfcs(lensfc), alffcs(lensfc,2)
321  REAL :: cvfcs(lensfc), cvtfcs(lensfc)
322  REAL :: cvbfcs(lensfc), tprcp(lensfc)
323  REAL :: srflag(lensfc), swdfcs(lensfc)
324  REAL :: slcfcs(lensfc,lsoil), vmxfcs(lensfc)
325  REAL :: vmnfcs(lensfc), t2m(lensfc)
326  REAL :: q2m(lensfc), slpfcs(lensfc)
327  REAL :: absfcs(lensfc), orog_uf(lensfc)
328  REAL :: ustar(lensfc)
329  REAL :: fmm(lensfc), fhh(lensfc)
330  REAL :: rla(lensfc), rlo(lensfc)
331  REAL(KIND=4) :: zsoil(lsoil)
332  REAL :: sig1t(lensfc)
335  REAL, ALLOCATABLE :: slifcs_fg(:)
336 
337  TYPE(nsst_data) :: nsst
338  real, dimension(idim,jdim) :: tf_clm,tf_trd,sal_clm
339  real, dimension(lensfc) :: tf_clm_tile,tf_trd_tile,sal_clm_tile
340 
341 !--------------------------------------------------------------------------------
342 ! GSI_FILE is the path/name of the gaussian GSI file which contains NSST
343 ! increments.
344 !--------------------------------------------------------------------------------
345 
346  DATA gsi_file/'NULL'/
347 
348  namelist/namsfcd/ gsi_file
349 
350  sig1t = 0.0 ! Not a dead start!
351 
352  input_nml_file = "NULL"
353 
354  CALL baopenr(37, "fort.37", ierr)
355  READ (37, nml=namsfcd)
356  WRITE(6,namsfcd)
357 
358  print*
359  print*,'IN ROUTINE SFCDRV,IDIM=',idim,'JDIM=',jdim,'FH=',fh
360 
361 !--------------------------------------------------------------------------------
362 ! READ THE OROGRAPHY AND GRID POINT LAT/LONS FOR THE CUBED-SPHERE TILE.
363 !--------------------------------------------------------------------------------
364 
365  CALL read_lat_lon_orog(rla,rlo,orog,orog_uf,tile_num,idim,jdim,lensfc)
366 
367  DO i = 1, idim
368  idum(i,:) = i
369  ENDDO
370 
371  i_index = reshape(idum, (/lensfc/))
372 
373  DO i = 1, jdim
374  idum(:,i) = i
375  ENDDO
376 
377  j_index = reshape(idum, (/lensfc/) )
378 
379  IF (do_nsst) THEN
380  print*
381  print*,"WILL PROCESS NSST RECORDS."
382  ALLOCATE(nsst%C_0(lensfc))
383  ALLOCATE(nsst%C_D(lensfc))
384  ALLOCATE(nsst%D_CONV(lensfc))
385  ALLOCATE(nsst%DT_COOL(lensfc))
386  ALLOCATE(nsst%IFD(lensfc))
387  ALLOCATE(nsst%QRAIN(lensfc))
388  ALLOCATE(nsst%TREF(lensfc))
389  ALLOCATE(nsst%TFINC(lensfc))
390  ALLOCATE(nsst%W_0(lensfc))
391  ALLOCATE(nsst%W_D(lensfc))
392  ALLOCATE(nsst%XS(lensfc))
393  ALLOCATE(nsst%XT(lensfc))
394  ALLOCATE(nsst%XTTS(lensfc))
395  ALLOCATE(nsst%XU(lensfc))
396  ALLOCATE(nsst%XV(lensfc))
397  ALLOCATE(nsst%XZ(lensfc))
398  ALLOCATE(nsst%XZTS(lensfc))
399  ALLOCATE(nsst%Z_C(lensfc))
400  ALLOCATE(nsst%ZM(lensfc))
401  ALLOCATE(slifcs_fg(lensfc))
402  ENDIF
403 
404 !--------------------------------------------------------------------------------
405 ! READ THE INPUT SURFACE DATA ON THE CUBED-SPHERE TILE.
406 !--------------------------------------------------------------------------------
407 
408  CALL read_data(tsffcs,smcfcs,snofcs,stcfcs,tg3fcs,zorfcs, &
409  cvfcs,cvbfcs,cvtfcs,albfcs,slifcs, &
410  vegfcs,cnpfcs,f10m,vetfcs,sotfcs, &
411  alffcs,ustar,fmm,fhh,sihfcs,sicfcs, &
412  sitfcs,tprcp,srflag,swdfcs,vmnfcs, &
413  vmxfcs,slcfcs,slpfcs,absfcs,t2m,q2m, &
414  slmask,zsoil,lsoil,lensfc,do_nsst,nsst)
415 
416  IF (use_ufo) THEN
417  print*
418  print*,'USE UNFILTERED OROGRAPHY.'
419  ELSE
420  orog_uf = 0.0
421  ENDIF
422 
423  DO i=1,lensfc
424  aisfcs(i) = 0.
425  IF(nint(slifcs(i)).EQ.2) aisfcs(i) = 1.
426  ENDDO
427 
428  IF (do_nsst) THEN
429  IF (adjt_nst_only) THEN
430  print*
431  print*,"FIRST GUESS MASK ADJUSTED BY IFD RECORD"
432  slifcs_fg = slifcs
433  WHERE(nint(nsst%IFD) == 3) slifcs_fg = 2.0
434  ELSE
435  print*
436  print*,"SAVE FIRST GUESS MASK"
437  slifcs_fg = slifcs
438  ENDIF
439  ENDIF
440 
441 !--------------------------------------------------------------------------------
442 ! UPDATE SURFACE FIELDS.
443 !--------------------------------------------------------------------------------
444 
445  IF (.NOT. adjt_nst_only) THEN
446  print*
447  print*,"CALL SFCCYCLE TO UPDATE SURFACE FIELDS."
448  CALL sfccycle(lugb,lensfc,lsoil,sig1t,deltsfc, &
449  iy,im,id,ih,fh,rla,rlo, &
450  slmask,orog, orog_uf, use_ufo, do_nsst, &
451  sihfcs,sicfcs,sitfcs,swdfcs,slcfcs, &
452  vmnfcs,vmxfcs,slpfcs,absfcs, &
453  tsffcs,snofcs,zorfcs,albfcs,tg3fcs, &
454  cnpfcs,smcfcs,stcfcs,slifcs,aisfcs, &
455  vegfcs,vetfcs,sotfcs,alffcs, &
456  cvfcs,cvbfcs,cvtfcs,myrank,nlunit, &
457  sz_nml, input_nml_file, &
458  ialb,isot,ivegsrc,tile_num,i_index,j_index)
459  ENDIF
460 
461 !--------------------------------------------------------------------------------
462 ! IF RUNNING WITH NSST, READ IN GSI FILE WITH THE UPDATED INCREMENTS (ON THE
463 ! GAUSSIAN GRID), INTERPOLATE INCREMENTS TO THE CUBED-SPHERE TILE, AND PERFORM
464 ! REQUIRED ADJUSTMENTS AND QC.
465 !--------------------------------------------------------------------------------
466 
467  IF (do_nsst) THEN
468  IF (gsi_file == "NULL") THEN
469  print*
470  print*,"NO GSI FILE. ADJUST IFD FOR FORMER ICE POINTS."
471  DO i = 1, lensfc
472  IF (nint(slifcs_fg(i)) == 2 .AND. nint(slifcs(i)) == 0) THEN
473  nsst%IFD(i) = 3.0
474  ENDIF
475  ENDDO
476  nsst%TFINC = 0.0
477  ELSE
478  print*
479  print*,"ADJUST TREF FROM GSI INCREMENT"
480 !
481 ! Get tf climatology at the time
482 !
483  call get_tf_clm(rla,rlo,jdim,idim,iy,im,id,ih,tf_clm,tf_trd)
484  tf_clm_tile(:) = reshape(tf_clm, (/lensfc/) )
485  tf_trd_tile(:) = reshape(tf_trd, (/lensfc/) )
486 !
487 ! Get salinity climatology at the time
488 !
489  call get_sal_clm(rla,rlo,jdim,idim,iy,im,id,ih,sal_clm)
490  sal_clm_tile(:) = reshape(sal_clm, (/lensfc/) )
491 !
492 ! read tf analysis increment generated by GSI
493 !
494  CALL read_gsi_data(gsi_file)
495 !
496 ! update foundation & surface temperature for NSST
497 !
498  CALL adjust_nsst(rla,rlo,slifcs,slifcs_fg,tsffcs,sitfcs,sicfcs,stcfcs, &
499  nsst,lensfc,lsoil,idim,jdim,zsea1,zsea2,im,id,deltsfc, &
500  tf_clm_tile,tf_trd_tile,sal_clm_tile)
501  ENDIF
502  ENDIF
503 
504 !--------------------------------------------------------------------------------
505 ! WRITE OUT UPDATED SURFACE DATA ON THE CUBED-SPHERE TILE.
506 !--------------------------------------------------------------------------------
507 
508  CALL write_data(slifcs,tsffcs,snofcs,tg3fcs,zorfcs, &
509  albfcs,alffcs,vegfcs,cnpfcs,f10m, &
510  t2m,q2m,vetfcs,sotfcs,ustar,fmm,fhh, &
511  sicfcs,sihfcs,sitfcs, &
512  tprcp,srflag,swdfcs, &
513  vmnfcs,vmxfcs,slpfcs,absfcs, &
514  slcfcs,smcfcs,stcfcs, &
515  idim,jdim,lensfc,lsoil,do_nsst,nsst)
516 
517  IF (do_nsst) THEN
518  DEALLOCATE(nsst%C_0)
519  DEALLOCATE(nsst%C_D)
520  DEALLOCATE(nsst%D_CONV)
521  DEALLOCATE(nsst%DT_COOL)
522  DEALLOCATE(nsst%IFD)
523  DEALLOCATE(nsst%QRAIN)
524  DEALLOCATE(nsst%TREF)
525  DEALLOCATE(nsst%TFINC)
526  DEALLOCATE(nsst%W_0)
527  DEALLOCATE(nsst%W_D)
528  DEALLOCATE(nsst%XS)
529  DEALLOCATE(nsst%XT)
530  DEALLOCATE(nsst%XTTS)
531  DEALLOCATE(nsst%XU)
532  DEALLOCATE(nsst%XV)
533  DEALLOCATE(nsst%XZ)
534  DEALLOCATE(nsst%XZTS)
535  DEALLOCATE(nsst%Z_C)
536  DEALLOCATE(nsst%ZM)
537  DEALLOCATE(slifcs_fg)
538  ENDIF
539 
540  RETURN
541 
542  END SUBROUTINE sfcdrv
543 
575  SUBROUTINE adjust_nsst(RLA,RLO,SLMSK_TILE,SLMSK_FG_TILE,SKINT_TILE,&
576  sicet_tile,sice_tile,soilt_tile,nsst,lensfc,lsoil, &
577  idim,jdim,zsea1,zsea2,mon,day,deltsfc, &
578  tf_clm_tile,tf_trd_tile,sal_clm_tile)
579 
580  USE gdswzd_mod
581  USE read_write_data, ONLY : idim_gaus, jdim_gaus, &
582  slmsk_gaus, dtref_gaus, &
583  nsst_data
584 
585  USE mpi
586 
587  IMPLICIT NONE
588 
589  INTEGER, INTENT(IN) :: lensfc, lsoil, idim, jdim, mon, day
590 
591  REAL, INTENT(IN) :: slmsk_tile(lensfc), slmsk_fg_tile(lensfc)
592  real, intent(in) :: tf_clm_tile(lensfc),tf_trd_tile(lensfc),sal_clm_tile(lensfc)
593  REAL, INTENT(IN) :: zsea1, zsea2, deltsfc
594  REAL, INTENT(INOUT) :: rla(lensfc), rlo(lensfc), skint_tile(lensfc)
595  REAL, INTENT(INOUT) :: sicet_tile(lensfc),sice_tile(lensfc),soilt_tile(lensfc,lsoil)
596 
597  TYPE(nsst_data) :: nsst
598 
599  REAL, PARAMETER :: tmax=313.0,tzero=273.16
600 
601  INTEGER :: iopt, nret, kgds_gaus(200)
602  INTEGER :: igaus, jgaus, ij, ii, jj, iii, jjj, krad
603  INTEGER :: istart, iend, jstart, jend
604  INTEGER :: mask_tile, mask_fg_tile
605  INTEGER :: itile, jtile
606  INTEGER :: max_search, j, ierr
607  INTEGER :: igausp1, jgausp1
608  integer :: nintp,nsearched,nice,nland
609  integer :: nfill,nfill_tice,nfill_clm
610  integer :: nset_thaw,nset_thaw_s,nset_thaw_i,nset_thaw_c
611 
612  INTEGER, ALLOCATABLE :: id1(:,:), id2(:,:), jdc(:,:)
613 
614  LOGICAL :: is_ice
615 
616  real :: tfreez
617  REAL :: tref_save,wsum,tf_ice,tf_thaw
618  REAL :: fill, dtzm, gaus_res_km, dtref
619  REAL, ALLOCATABLE :: xpts(:), ypts(:), lats(:), lons(:)
620  REAL, ALLOCATABLE :: dum2d(:,:), lats_rad(:), lons_rad(:)
621  REAL, ALLOCATABLE :: agrid(:,:,:), s2c(:,:,:)
622 
623  kgds_gaus = 0
624  kgds_gaus(1) = 4 ! OCT 6 - TYPE OF GRID (GAUSSIAN)
625  kgds_gaus(2) = idim_gaus ! OCT 7-8 - # PTS ON LATITUDE CIRCLE
626  kgds_gaus(3) = jdim_gaus
627  kgds_gaus(4) = 90000 ! OCT 11-13 - LAT OF ORIGIN
628  kgds_gaus(5) = 0 ! OCT 14-16 - LON OF ORIGIN
629  kgds_gaus(6) = 128 ! OCT 17 - RESOLUTION FLAG
630  kgds_gaus(7) = -90000 ! OCT 18-20 - LAT OF EXTREME POINT
631  kgds_gaus(8) = nint(-360000./float(idim_gaus)) ! OCT 21-23 - LON OF EXTREME POINT
632  kgds_gaus(9) = nint((360.0 / float(idim_gaus))*1000.0)
633  ! OCT 24-25 - LONGITUDE DIRECTION INCR.
634  kgds_gaus(10) = jdim_gaus/2 ! OCT 26-27 - NUMBER OF CIRCLES POLE TO EQUATOR
635  kgds_gaus(12) = 255 ! OCT 29 - RESERVED
636  kgds_gaus(20) = 255 ! OCT 5 - NOT USED, SET TO 255
637 
638  print*
639  print*,'ADJUST NSST USING GSI INCREMENTS ON GAUSSIAN GRID'
640 
641 !----------------------------------------------------------------------
642 ! CALL GDSWZD TO COMPUTE THE LAT/LON OF EACH GSI GAUSSIAN GRID POINT.
643 !----------------------------------------------------------------------
644 
645  iopt = 0
646  fill = -9999.
647  ALLOCATE(xpts(idim_gaus*jdim_gaus))
648  ALLOCATE(ypts(idim_gaus*jdim_gaus))
649  ALLOCATE(lats(idim_gaus*jdim_gaus))
650  ALLOCATE(lons(idim_gaus*jdim_gaus))
651  xpts = fill
652  ypts = fill
653  lats = fill
654  lons = fill
655 
656  CALL gdswzd(kgds_gaus,iopt,(idim_gaus*jdim_gaus),fill,xpts,ypts,lons,lats,nret)
657 
658  IF (nret /= (idim_gaus*jdim_gaus)) THEN
659  print*,'FATAL ERROR: PROBLEM IN GDSWZD. STOP.'
660  CALL mpi_abort(mpi_comm_world, 12, ierr)
661  ENDIF
662 
663  DEALLOCATE (xpts, ypts)
664 
665  ALLOCATE(dum2d(idim_gaus,jdim_gaus))
666  dum2d = reshape(lats, (/idim_gaus,jdim_gaus/) )
667  DEALLOCATE(lats)
668 
669  ALLOCATE(lats_rad(jdim_gaus))
670  DO j = 1, jdim_gaus
671  lats_rad(j) = dum2d(1,jdim_gaus-j+1) * 3.1415926 / 180.0
672  ENDDO
673 
674  dum2d = reshape(lons, (/idim_gaus,jdim_gaus/) )
675  DEALLOCATE(lons)
676  ALLOCATE(lons_rad(idim_gaus))
677  lons_rad = dum2d(:,1) * 3.1415926 / 180.0
678 
679  DEALLOCATE(dum2d)
680 
681  ALLOCATE(agrid(idim,jdim,2))
682  agrid(:,:,1) = reshape(rlo, (/idim,jdim/) )
683  agrid(:,:,2) = reshape(rla, (/idim,jdim/) )
684  agrid = agrid * 3.1415926 / 180.0
685 
686  ALLOCATE(id1(idim,jdim))
687  ALLOCATE(id2(idim,jdim))
688  ALLOCATE(jdc(idim,jdim))
689  ALLOCATE(s2c(idim,jdim,4))
690 
691 !----------------------------------------------------------------------
692 ! COMPUTE BILINEAR WEIGHTS FOR EACH MODEL POINT FROM THE NEAREST
693 ! FOUR GSI/GAUSSIAN POINTS. DOES NOT ACCOUNT FOR MASK. THAT
694 ! HAPPENS LATER.
695 !----------------------------------------------------------------------
696 
697  CALL remap_coef( 1, idim, 1, jdim, idim_gaus, jdim_gaus, &
698  lons_rad, lats_rad, id1, id2, jdc, s2c, agrid )
699 
700  DEALLOCATE(lons_rad, lats_rad, agrid)
701 
702 !----------------------------------------------------------------------
703 ! THE MAXIMUM DISTANCE TO SEARCH IS 500 KM. HOW MANY GAUSSIAN
704 ! GRID LENGTHS IS THAT?
705 !----------------------------------------------------------------------
706 
707  gaus_res_km = 360.0 / idim_gaus * 111.0
708  max_search = ceiling(500.0/gaus_res_km)
709 
710  print*
711  print*,'MAXIMUM SEARCH IS ',max_search, ' GAUSSIAN POINTS.'
712  print*
713 
714 !
715 ! Initialize variables for counts statitics to be zeros
716 !
717  nintp = 0
718  nset_thaw = 0
719  nset_thaw_s = 0
720  nset_thaw_i = 0
721  nset_thaw_c = 0
722  nsearched = 0
723  nfill = 0
724  nfill_tice = 0
725  nfill_clm = 0
726  nice = 0
727  nland = 0
728 !----------------------------------------------------------------------
729 ! TREF INCREMENT WILL BE OUTPUT. INITIALIZE TO ZERO.
730 !----------------------------------------------------------------------
731 
732  nsst%TFINC = 0.0
733 
734  ij_loop : DO ij = 1, lensfc
735 
736  mask_tile = nint(slmsk_tile(ij))
737  mask_fg_tile = nint(slmsk_fg_tile(ij))
738 
739 !
740 ! when sea ice exists, get salinity dependent water temperature
741 !
742  tf_ice = tfreez(sal_clm_tile(ij)) + tzero
743 !----------------------------------------------------------------------
744 ! SKIP LAND POINTS. NSST NOT APPLIED AT LAND.
745 !----------------------------------------------------------------------
746 
747  IF (mask_tile == 1) THEN
748  nland = nland + 1
749  cycle ij_loop
750  ENDIF
751 
752 !
753 ! these are ice points. set tref to tf_ice and update tmpsfc.
754 !
755  if (mask_tile == 2) then
756  nsst%tref(ij)=tf_ice ! water part tmp set
757  skint_tile(ij)=(1.0-sice_tile(ij))*nsst%tref(ij)+sice_tile(ij)*sicet_tile(ij)
758  nice = nice + 1
759  cycle ij_loop
760  endif
761 
762 !
763 ! Get i,j index on array of (idim,jdim) from known ij
764 !
765  jtile = (ij-1) / idim + 1
766  itile = mod(ij,idim)
767  IF (itile==0) itile = idim
768 
769 !----------------------------------------------------------------------
770 ! IF THE MODEL POINT WAS ICE COVERED, BUT IS NOW OPEN WATER, SET
771 ! TREF TO searched adjascent open water onea, if failed the search, set to
772 ! weighted average of tf_ice and tf_clm. For NSST vars, set xz TO '30' AND ALL OTHER FIELDS TO ZERO.
773 !----------------------------------------------------------------------
774 
775  IF (mask_fg_tile == 2 .AND. mask_tile == 0) THEN
776 !
777 ! set background for the thaw (just melted water) situation
778 !
779  call tf_thaw_set(nsst%tref,nint(slmsk_fg_tile),itile,jtile,tf_ice,tf_clm_tile(ij),tf_thaw,idim,jdim, &
780  nset_thaw_s,nset_thaw_i,nset_thaw_c)
781  call nsst_water_reset(nsst,ij,tf_thaw)
782  nset_thaw = nset_thaw + 1
783  ENDIF
784 
785 !----------------------------------------------------------------------
786 ! THESE ARE POINTS THAT ARE OPEN WATER AND WERE OPEN WATER PRIOR
787 ! TO ANY ICE UPDATE BY SFCCYCLE. UPDATE TREF AND SKIN TEMP.
788 ! AT OPEN WATER POINTS, THE SEA ICE TEMPERATURE (SICET_TILE) AND
789 ! SOIL COLUMN TEMPERATURE (SOILT_TILE) ARE SET TO THE SKIN TEMP.
790 ! IT IS SIMPLY A FILLER VALUE. THESE FIELDS ARE NOT USED AT
791 ! OPEN WATER POINTS.
792 !----------------------------------------------------------------------
793 !----------------------------------------------------------------------
794 ! SEE IF ANY OF THE NEAREST GSI POINTS MASK AREA OPEN WATER.
795 ! IF SO, APPLY NSST INCREMENT USING BILINEAR INTERPOLATION.
796 !----------------------------------------------------------------------
797 
798  igaus = id1(itile,jtile)
799  jgaus = jdc(itile,jtile)
800  igausp1 = id2(itile,jtile)
801  jgausp1 = jdc(itile,jtile)+1
802 
803  IF (slmsk_gaus(igaus,jgaus) == 0 .OR. &
804  slmsk_gaus(igausp1,jgaus) == 0 .OR. &
805  slmsk_gaus(igausp1,jgausp1) == 0 .OR. &
806  slmsk_gaus(igaus,jgausp1) == 0) THEN
807 
808  dtref = 0.0
809  wsum = 0.0
810 
811  IF (slmsk_gaus(igaus,jgaus) == 0) THEN
812  dtref = dtref + (s2c(itile,jtile,1) * dtref_gaus(igaus,jgaus))
813  wsum = wsum + s2c(itile,jtile,1)
814  ENDIF
815 
816  IF (slmsk_gaus(igausp1,jgaus) == 0) THEN
817  dtref = dtref + (s2c(itile,jtile,2) * dtref_gaus(igausp1,jgaus))
818  wsum = wsum + s2c(itile,jtile,2)
819  ENDIF
820 
821  IF (slmsk_gaus(igausp1,jgausp1) == 0) THEN
822  dtref = dtref + (s2c(itile,jtile,3) * dtref_gaus(igausp1,jgausp1))
823  wsum = wsum + s2c(itile,jtile,3)
824  ENDIF
825 
826  IF (slmsk_gaus(igaus,jgausp1) == 0) THEN
827  dtref = dtref + (s2c(itile,jtile,4) * dtref_gaus(igaus,jgausp1))
828  wsum = wsum + s2c(itile,jtile,4)
829  ENDIF
830 
831  nintp = nintp + 1
832  dtref = dtref / wsum
833 
834  tref_save = nsst%TREF(ij)
835  nsst%TREF(ij) = nsst%TREF(ij) + dtref
836  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
837  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
838  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
839 
840  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
841  nsst%Z_C(ij),zsea1,zsea2,dtzm)
842 
843  skint_tile(ij) = nsst%TREF(ij) + dtzm
844  skint_tile(ij) = max(skint_tile(ij), tf_ice)
845  skint_tile(ij) = min(skint_tile(ij), tmax)
846 
847  sicet_tile(ij) = skint_tile(ij)
848  soilt_tile(ij,:) = skint_tile(ij)
849 
850 !----------------------------------------------------------------------
851 ! NO NEARBY GSI/GAUSSIAN OPEN WATER POINTS. PERFORM A SPIRAL SEARCH TO
852 ! FIND NEAREST NON-LAND POINT ON GSI/GAUSSIAN GRID.
853 !----------------------------------------------------------------------
854 
855  ELSE
856 
857  is_ice = .false.
858 
859  DO krad = 1, max_search
860 
861  istart = igaus - krad
862  iend = igaus + krad
863  jstart = jgaus - krad
864  jend = jgaus + krad
865 
866  DO jj = jstart, jend
867  DO ii = istart, iend
868 
869  IF((jj == jstart) .OR. (jj == jend) .OR. &
870  (ii == istart) .OR. (ii == iend)) THEN
871 
872  IF ((jj >= 1) .AND. (jj <= jdim_gaus)) THEN
873 
874  jjj = jj
875  IF (ii <= 0) THEN
876  iii = idim_gaus + ii
877  ELSE IF (ii >= (idim_gaus+1)) THEN
878  iii = ii - idim_gaus
879  ELSE
880  iii = ii
881  END IF
882 
883 !----------------------------------------------------------------------
884 ! SEE IF NEARBY POINTS ARE SEA ICE. IF THEY ARE, AND THE SEARCH FOR
885 ! A GAUSSIAN GRID OPEN WATER POINT FAILS, THEN TREF WILL BE SET TO
886 ! FREEZING BELOW.
887 !----------------------------------------------------------------------
888 
889  IF (krad <= 2 .AND. slmsk_gaus(iii,jjj) == 2) is_ice = .true.
890 
891  IF (slmsk_gaus(iii,jjj) == 0) THEN
892 
893 ! PRINT*,'MISMATCH AT TILE POINT ',ITILE,JTILE
894 ! PRINT*,'UPDATE TREF USING GSI INCREMENT AT ',III,JJJ,DTREF_GAUS(III,JJJ)
895  nsearched = nsearched + 1
896 
897  tref_save = nsst%TREF(ij)
898  nsst%TREF(ij ) = nsst%TREF(ij) + dtref_gaus(iii,jjj)
899  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
900  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
901  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
902 
903  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
904  nsst%Z_C(ij),zsea1,zsea2,dtzm)
905 
906  skint_tile(ij) = nsst%TREF(ij) + dtzm
907  skint_tile(ij) = max(skint_tile(ij), tf_ice)
908  skint_tile(ij) = min(skint_tile(ij), tmax)
909 
910  sicet_tile(ij) = skint_tile(ij)
911  soilt_tile(ij,:) = skint_tile(ij)
912  cycle ij_loop
913 
914  ENDIF ! GSI/Gaussian mask is open water
915 
916  ENDIF
917 
918  ENDIF
919 
920  ENDDO
921  ENDDO
922 
923  ENDDO ! KRAD LOOP
924 
925 !----------------------------------------------------------------------
926 ! THE SEARCH FAILED. IF THERE IS NEARBY ICE, SET TREF TO FREEZING.
927 ! ELSE UPDATE TREF BASED ON THE ANNUAL SST CYCLE.
928 !----------------------------------------------------------------------
929 
930 ! PRINT*,'WARNING !!!!!! SEARCH FAILED AT TILE POINT ',ITILE,JTILE
931 
932  nfill = nfill + 1
933  IF (is_ice) THEN
934  nsst%TREF(ij) = tf_ice
935 ! PRINT*,"NEARBY ICE. SET TREF TO FREEZING"
936  nfill_tice = nfill_tice + 1
937  ELSE
938  tref_save = nsst%TREF(ij)
939  nsst%TREF(ij) = nsst%TREF(ij) + tf_trd_tile(ij)
940  nsst%TREF(ij) = max(nsst%TREF(ij), tf_ice)
941  nsst%TREF(ij) = min(nsst%TREF(ij), tmax)
942  nsst%TFINC(ij) = nsst%TREF(ij) - tref_save
943 ! PRINT*,'UPDATE TREF FROM SST CLIMO ',DTREF
944  nfill_clm = nfill_clm + 1
945  ENDIF
946 
947  CALL dtzm_point(nsst%XT(ij),nsst%XZ(ij),nsst%DT_COOL(ij), &
948  nsst%Z_C(ij),zsea1,zsea2,dtzm)
949 
950  skint_tile(ij) = nsst%TREF(ij) + dtzm
951  skint_tile(ij) = max(skint_tile(ij), tf_ice)
952  skint_tile(ij) = min(skint_tile(ij), tmax)
953 
954  sicet_tile(ij) = skint_tile(ij)
955  soilt_tile(ij,:) = skint_tile(ij)
956 
957  ENDIF ! NEARBY GAUSSIAN POINTS ARE OPEN WATER?
958 
959  ENDDO ij_loop
960 
961  write(*,'(a)') 'statistics of grids number processed for tile : '
962  write(*,'(a,I8)') ' nintp = ',nintp
963  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
964  write(*,'(a,I8)') ' nsearched = ',nsearched
965  write(*,'(a,3I6)') ' nfill,nfill_tice,nfill_clm = ',nfill,nfill_tice,nfill_clm
966  write(*,'(a,I8)') ' nice = ',nice
967  write(*,'(a,I8)') ' nland = ',nland
968 
969  DEALLOCATE(id1, id2, jdc, s2c)
970 
971  END SUBROUTINE adjust_nsst
972 
983  SUBROUTINE climo_trend(LATITUDE, MON, DAY, DELTSFC, DTREF)
984  IMPLICIT NONE
985 
986  INTEGER, INTENT(IN) :: mon, day
987 
988  REAL, INTENT(IN) :: latitude, deltsfc
989  REAL, INTENT(OUT) :: dtref
990 
991  INTEGER :: num_days(12), mon2, mon1
992 
993  REAL, TARGET :: sst_80_90(12)
994  REAL, TARGET :: sst_70_80(12)
995  REAL, TARGET :: sst_60_70(12)
996  REAL, TARGET :: sst_50_60(12)
997  REAL, TARGET :: sst_40_50(12)
998  REAL, TARGET :: sst_30_40(12)
999  REAL, TARGET :: sst_20_30(12)
1000  REAL, TARGET :: sst_10_20(12)
1001  REAL, TARGET :: sst_00_10(12)
1002  REAL, TARGET :: sst_m10_00(12)
1003  REAL, TARGET :: sst_m20_m10(12)
1004  REAL, TARGET :: sst_m30_m20(12)
1005  REAL, TARGET :: sst_m40_m30(12)
1006  REAL, TARGET :: sst_m50_m40(12)
1007  REAL, TARGET :: sst_m60_m50(12)
1008  REAL, TARGET :: sst_m70_m60(12)
1009  REAL, TARGET :: sst_m80_m70(12)
1010  REAL, TARGET :: sst_m90_m80(12)
1011 
1012  REAL, POINTER :: sst(:)
1013 
1014  DATA num_days /31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31/
1015 
1016  DATA sst_80_90 /271.466, 271.458, 271.448, 271.445, 271.519, 271.636, &
1017  272.023, 272.066, 272.001, 271.698, 271.510, 271.472/
1018 
1019  DATA sst_70_80 /272.149, 272.103, 272.095, 272.126, 272.360, 272.988, &
1020  274.061, 274.868, 274.415, 273.201, 272.468, 272.268/
1021 
1022  DATA sst_60_70 /274.240, 274.019, 273.988, 274.185, 275.104, 276.875, &
1023  279.005, 280.172, 279.396, 277.586, 275.818, 274.803/
1024 
1025  DATA sst_50_60 /277.277, 276.935, 277.021, 277.531, 279.100, 281.357, &
1026  283.735, 285.171, 284.399, 282.328, 279.918, 278.199/
1027 
1028  DATA sst_40_50 /281.321, 280.721, 280.850, 281.820, 283.958, 286.588, &
1029  289.195, 290.873, 290.014, 287.652, 284.898, 282.735/
1030 
1031  DATA sst_30_40 /289.189, 288.519, 288.687, 289.648, 291.547, 293.904, &
1032  296.110, 297.319, 296.816, 295.225, 292.908, 290.743/
1033 
1034  DATA sst_20_30 /294.807, 294.348, 294.710, 295.714, 297.224, 298.703, &
1035  299.682, 300.127, 300.099, 299.455, 297.953, 296.177/
1036 
1037  DATA sst_10_20 /298.878, 298.720, 299.033, 299.707, 300.431, 300.709, &
1038  300.814, 300.976, 301.174, 301.145, 300.587, 299.694/
1039 
1040  DATA sst_00_10 /300.415, 300.548, 300.939, 301.365, 301.505, 301.141, &
1041  300.779, 300.660, 300.818, 300.994, 300.941, 300.675/
1042 
1043  DATA sst_m10_00 /300.226, 300.558, 300.914, 301.047, 300.645, 299.870, &
1044  299.114, 298.751, 298.875, 299.294, 299.721, 299.989/
1045 
1046  DATA sst_m20_m10 /299.547, 299.985, 300.056, 299.676, 298.841, 297.788, &
1047  296.893, 296.491, 296.687, 297.355, 298.220, 298.964/
1048 
1049  DATA sst_m30_m20 /297.524, 298.073, 297.897, 297.088, 295.846, 294.520, &
1050  293.525, 293.087, 293.217, 293.951, 295.047, 296.363/
1051 
1052  DATA sst_m40_m30 /293.054, 293.765, 293.468, 292.447, 291.128, 289.781, &
1053  288.773, 288.239, 288.203, 288.794, 289.947, 291.553/
1054 
1055  DATA sst_m50_m40 /285.052, 285.599, 285.426, 284.681, 283.761, 282.826, &
1056  282.138, 281.730, 281.659, 281.965, 282.768, 283.961/
1057 
1058  DATA sst_m60_m50 /277.818, 278.174, 277.991, 277.455, 276.824, 276.229, &
1059  275.817, 275.585, 275.560, 275.687, 276.142, 276.968/
1060 
1061  DATA sst_m70_m60 /273.436, 273.793, 273.451, 272.813, 272.349, 272.048, &
1062  271.901, 271.838, 271.845, 271.889, 272.080, 272.607/
1063 
1064  DATA sst_m80_m70 /271.579, 271.578, 271.471, 271.407, 271.392, 271.391, &
1065  271.390, 271.391, 271.394, 271.401, 271.422, 271.486/
1066 
1067  DATA sst_m90_m80 /271.350, 271.350, 271.350, 271.350, 271.350, 271.350, &
1068  271.350, 271.350, 271.350, 271.350, 271.350, 271.350/
1069 
1070  nullify(sst)
1071  IF (latitude > 80.0) THEN
1072  sst => sst_80_90
1073  ELSEIF (latitude > 70.0) THEN
1074  sst => sst_70_80
1075  ELSEIF (latitude > 60.0) THEN
1076  sst => sst_60_70
1077  ELSEIF (latitude > 50.0) THEN
1078  sst => sst_50_60
1079  ELSEIF (latitude > 40.0) THEN
1080  sst => sst_40_50
1081  ELSEIF (latitude > 30.0) THEN
1082  sst => sst_30_40
1083  ELSEIF (latitude > 20.0) THEN
1084  sst => sst_20_30
1085  ELSEIF (latitude > 10.0) THEN
1086  sst => sst_10_20
1087  ELSEIF (latitude > 0.0) THEN
1088  sst => sst_00_10
1089  ELSEIF (latitude > -10.0) THEN
1090  sst => sst_m10_00
1091  ELSEIF (latitude > -20.0) THEN
1092  sst => sst_m20_m10
1093  ELSEIF (latitude > -30.0) THEN
1094  sst => sst_m30_m20
1095  ELSEIF (latitude > -40.0) THEN
1096  sst => sst_m40_m30
1097  ELSEIF (latitude > -50.0) THEN
1098  sst => sst_m50_m40
1099  ELSEIF (latitude > -60.0) THEN
1100  sst => sst_m60_m50
1101  ELSEIF (latitude > -70.0) THEN
1102  sst => sst_m70_m60
1103  ELSEIF (latitude > -80.0) THEN
1104  sst => sst_m80_m70
1105  ELSE
1106  sst => sst_m90_m80
1107  END IF
1108 
1109  IF (day >= 15) THEN
1110  mon2 = mon+1
1111  IF(mon2 == 13) mon2 = 1
1112  mon1 = mon
1113  dtref = (sst(mon2) - sst(mon1)) / num_days(mon1)
1114  ELSE
1115  mon1 = mon - 1
1116  IF (mon1 == 0) mon1=12
1117  mon2 = mon
1118  dtref = (sst(mon2) - sst(mon1)) / num_days(mon1)
1119  ENDIF
1120 
1121  dtref = dtref * (deltsfc / 24.0)
1122 
1123  END SUBROUTINE climo_trend
1124 
1136  SUBROUTINE dtzm_point(XT,XZ,DT_COOL,ZC,Z1,Z2,DTZM)
1137  implicit none
1138 
1139  real, intent(in) :: xt,xz,dt_cool,zc,z1,z2
1140  real, intent(out) :: dtzm
1141 
1142  real, parameter :: zero = 0.0
1143  real, parameter :: one = 1.0
1144  real, parameter :: half = 0.5
1145  real :: dt_warm,dtw,dtc
1146 !
1147 ! get the mean warming in the range of z=z1 to z=z2
1148 !
1149  dtw = zero
1150  if ( xt > zero ) then
1151  dt_warm = (xt+xt)/xz ! Tw(0)
1152  if ( z1 < z2) then
1153  if ( z2 < xz ) then
1154  dtw = dt_warm*(one-(z1+z2)/(xz+xz))
1155  elseif ( z1 < xz .and. z2 >= xz ) then
1156  dtw = half*(one-z1/xz)*dt_warm*(xz-z1)/(z2-z1)
1157  endif
1158  elseif ( z1 == z2 ) then
1159  if ( z1 < xz ) then
1160  dtw = dt_warm*(one-z1/xz)
1161  endif
1162  endif
1163  endif
1164 !
1165 ! get the mean cooling in the range of z=z1 to z=z2
1166 !
1167  dtc = zero
1168  if ( zc > zero ) then
1169  if ( z1 < z2) then
1170  if ( z2 < zc ) then
1171  dtc = dt_cool*(one-(z1+z2)/(zc+zc))
1172  elseif ( z1 < zc .and. z2 >= zc ) then
1173  dtc = half*(one-z1/zc)*dt_cool*(zc-z1)/(z2-z1)
1174  endif
1175  elseif ( z1 == z2 ) then
1176  if ( z1 < zc ) then
1177  dtc = dt_cool*(one-z1/zc)
1178  endif
1179  endif
1180  endif
1181 
1182 !
1183 ! get the mean T departure from Tf in the range of z=z1 to z=z2
1184 !
1185  dtzm = dtw - dtc
1186 
1187  END SUBROUTINE dtzm_point
1188 
1214  SUBROUTINE remap_coef( is, ie, js, je,&
1215  im, jm, lon, lat, id1, id2, jdc, s2c, agrid )
1216 
1217  implicit none
1218  integer, intent(in):: is, ie, js, je
1219  integer, intent(in):: im, jm
1220  real, intent(in):: lon(im), lat(jm)
1221  real, intent(out):: s2c(is:ie,js:je,4)
1222  integer, intent(out), dimension(is:ie,js:je):: id1, id2, jdc
1223  real, intent(in):: agrid(is:ie,js:je,2)
1224  ! local:
1225  real :: rdlon(im)
1226  real :: rdlat(jm)
1227  real:: a1, b1
1228  real, parameter :: pi = 3.1415926
1229  integer i,j, i1, i2, jc, i0, j0
1230  do i=1,im-1
1231  rdlon(i) = 1. / (lon(i+1) - lon(i))
1232  enddo
1233  rdlon(im) = 1. / (lon(1) + 2.*pi - lon(im))
1234 
1235  do j=1,jm-1
1236  rdlat(j) = 1. / (lat(j+1) - lat(j))
1237  enddo
1238 
1239  ! * Interpolate to cubed sphere cell center
1240  do 5000 j=js,je
1241 
1242  do i=is,ie
1243 
1244  if ( agrid(i,j,1)>lon(im) ) then
1245  i1 = im; i2 = 1
1246  a1 = (agrid(i,j,1)-lon(im)) * rdlon(im)
1247  elseif ( agrid(i,j,1)<lon(1) ) then
1248  i1 = im; i2 = 1
1249  a1 = (agrid(i,j,1)+2.*pi-lon(im)) * rdlon(im)
1250  else
1251  do i0=1,im-1
1252  if ( agrid(i,j,1)>=lon(i0) .and. agrid(i,j,1)<=lon(i0+1) ) then
1253  i1 = i0; i2 = i0+1
1254  a1 = (agrid(i,j,1)-lon(i1)) * rdlon(i0)
1255  go to 111
1256  endif
1257  enddo
1258  endif
1259 111 continue
1260 
1261  if ( agrid(i,j,2)<lat(1) ) then
1262  jc = 1
1263  b1 = 0.
1264  elseif ( agrid(i,j,2)>lat(jm) ) then
1265  jc = jm-1
1266  b1 = 1.
1267  else
1268  do j0=1,jm-1
1269  if ( agrid(i,j,2)>=lat(j0) .and. agrid(i,j,2)<=lat(j0+1) ) then
1270  jc = j0
1271  b1 = (agrid(i,j,2)-lat(jc)) * rdlat(jc)
1272  go to 222
1273  endif
1274  enddo
1275  endif
1276 222 continue
1277 
1278  if ( a1<0.0 .or. a1>1.0 .or. b1<0.0 .or. b1>1.0 ) then
1279  write(*,*) 'gid=', i,j,a1, b1
1280  endif
1281 
1282  s2c(i,j,1) = (1.-a1) * (1.-b1)
1283  s2c(i,j,2) = a1 * (1.-b1)
1284  s2c(i,j,3) = a1 * b1
1285  s2c(i,j,4) = (1.-a1) * b1
1286  id1(i,j) = i1
1287  id2(i,j) = i2
1288  jdc(i,j) = jc
1289  enddo !i-loop
1290 5000 continue ! j-loop
1291 
1292  END SUBROUTINE remap_coef
1293 
1313  subroutine tf_thaw_set(tf_ij,mask_ij,itile,jtile,tice,tclm,tf_thaw,nx,ny, &
1314  nset_thaw_s,nset_thaw_i,nset_thaw_c)
1315 
1316  real, dimension(nx*ny), intent(in) :: tf_ij
1317  integer, dimension(nx*ny), intent(in) :: mask_ij
1318  real, intent(in) :: tice,tclm
1319  integer, intent(in) :: itile,jtile,nx,ny
1320  real, intent(out) :: tf_thaw
1321  integer, intent(inout) :: nset_thaw_s,nset_thaw_i,nset_thaw_c
1322 ! Local
1323  real, parameter :: bmiss = -999.0
1324  real, dimension(nx,ny) :: tf
1325  integer, dimension(nx,ny) :: mask
1326  integer :: krad,max_search,istart,iend,jstart,jend
1327  integer :: ii,jj,iii,jjj
1328  logical :: is_ice
1329 
1330  max_search = 2
1331 
1332  mask(:,:) = reshape(mask_ij,(/nx,ny/) )
1333  tf(:,:) = reshape(tf_ij,(/nx,ny/) )
1334 
1335  tf_thaw = bmiss
1336 
1337  do krad = 1, max_search
1338 
1339  istart = itile - krad
1340  iend = itile + krad
1341  jstart = jtile - krad
1342  jend = jtile + krad
1343 
1344  do jj = jstart, jend
1345  do ii = istart, iend
1346 
1347 
1348  if ((jj == jstart) .or. (jj == jend) .or. &
1349  (ii == istart) .or. (ii == iend)) then
1350 
1351  if ((jj >= 1) .and. (jj <= ny)) then
1352  jjj = jj
1353  if (ii <= 0) then
1354  iii = nx + ii
1355  else if (ii >= (nx+1)) then
1356  iii = ii - nx
1357  else
1358  iii = ii
1359  endif
1360 
1361 !----------------------------------------------------------------------
1362 ! SEE IF NEARBY POINTS ARE SEA ICE. IF THEY ARE, AND THE SEARCH FOR
1363 ! A GAUSSIAN GRID OPEN WATER POINT FAILS, THEN TREF WILL BE SET TO
1364 ! FREEZING BELOW.
1365 !----------------------------------------------------------------------
1366  if (krad <= 2 .and. mask(iii,jjj) == 2) is_ice = .true.
1367 
1368  if (mask(iii,jjj) == 0) then
1369  tf_thaw = tf(iii,jjj)
1370  nset_thaw_s = nset_thaw_s + 1
1371  write(*,'(a,I4,2F9.3)') 'nset_thaw_s,tf(iii,jjj),tclm : ',nset_thaw_s,tf(iii,jjj),tclm
1372  go to 100
1373  endif ! tile mask is open water
1374 
1375  endif
1376  endif
1377  enddo
1378  enddo
1379  enddo ! krad loop
1380 
1381  100 continue
1382 
1383  if ( tf_thaw == bmiss ) then
1384  if (is_ice) then
1385  tf_thaw = tice
1386  nset_thaw_i = nset_thaw_i + 1
1387  write(*,'(a,I4,F9.3)') 'nset_thaw_i,tf_ice : ',nset_thaw_i,tice
1388  else
1389  tf_thaw = 0.8*tice+0.2*tclm
1390  nset_thaw_c = nset_thaw_c + 1
1391  write(*,'(a,I4,2F9.3)') 'nset_thaw_c,tf_ice,tclm : ',nset_thaw_c,tice,tclm
1392  endif
1393  endif
1394 
1395  end subroutine tf_thaw_set
1396 
1404  subroutine nsst_water_reset(nsst,ij,tf_thaw)
1405  use read_write_data, only : nsst_data
1406  implicit none
1407 
1408  integer, intent(in) :: ij
1409 
1410  real, intent(in) :: tf_thaw
1411 
1412  type(nsst_data), intent(inout) :: nsst
1413 
1414  nsst%c_0(ij) = 0.0
1415  nsst%c_d(ij) = 0.0
1416  nsst%d_conv(ij) = 0.0
1417  nsst%dt_cool(ij) = 0.0
1418  nsst%ifd(ij) = 0.0
1419  nsst%qrain(ij) = 0.0
1420  nsst%tref(ij) = tf_thaw
1421  nsst%w_0(ij) = 0.0
1422  nsst%w_d(ij) = 0.0
1423  nsst%xs(ij) = 0.0
1424  nsst%xt(ij) = 0.0
1425  nsst%xtts(ij) = 0.0
1426  nsst%xu(ij) = 0.0
1427  nsst%xv(ij) = 0.0
1428  nsst%xz(ij) = 30.0
1429  nsst%xzts(ij) = 0.0
1430  nsst%z_c(ij) = 0.0
1431  nsst%zm(ij) = 0.0
1432 
1433  end subroutine nsst_water_reset
1434 
1450 subroutine get_tf_clm(xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,tf_clm,tf_trd)
1451  use read_write_data, only : get_tf_clm_dim
1452 
1453  implicit none
1454 
1455  real, dimension(nx*ny), intent(in) :: xlats_ij
1456  real, dimension(nx*ny), intent(in) :: xlons_ij
1457  real, dimension(nx,ny), intent(out) :: tf_clm
1458  real, dimension(nx,ny), intent(out) :: tf_trd
1459  integer, intent(in) :: iy,im,id,ih,nx,ny
1460 ! local declare
1461  real, allocatable, dimension(:,:) :: tf_clm0 ! sst climatology at the valid time (nxc,nyc)
1462  real, allocatable, dimension(:,:) :: tf_trd0 ! 6-hourly sst climatology tendency valid at atime (nxc,nyc)
1463  real, allocatable, dimension(:) :: cxlats ! latitudes of sst climatology
1464  real, allocatable, dimension(:) :: cxlons ! longitudes of sst climatology
1465 
1466  real, dimension(nx*ny) :: tf_clm_ij ! sst climatology at target grids (nx*ny)
1467  real, dimension(nx*ny) :: tf_trd_ij ! 6-hourly sst climatology tendency
1468  real :: wei1,wei2
1469  integer :: nxc,nyc,mon1,mon2,i,j
1470  character (len=6), parameter :: fin_tf_clm='sstclm' ! sst climatology file name
1471 !
1472 ! get which two months used and their weights from atime
1473 !
1474  call get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1475 !
1476 ! get the dimensions of the sst climatology & allocate the related arrays
1477 !
1478  call get_tf_clm_dim(fin_tf_clm,nyc,nxc)
1479  allocate( tf_clm0(nxc,nyc),tf_trd0(nxc,nyc),cxlats(nyc),cxlons(nxc) )
1480 !
1481 ! get tf_clm at the analysis time from monthly climatology & cxlats, cxlons
1482 !
1483  call get_tf_clm_ta(tf_clm0,tf_trd0,cxlats,cxlons,nyc,nxc,mon1,mon2,wei1,wei2)
1484 !
1485 ! get tf_clm (nx by ny lat/lon) valid at atime
1486 !
1487  if ( nx == nxc .and. ny == nyc ) then
1488  tf_clm(:,:) = tf_clm0(:,:)
1489  tf_trd(:,:) = tf_trd0(:,:)
1490 ! write(*,'(a,2f9.3)') 'same dimensions, tf_clm, min : ',minval(tf_clm),maxval(tf_clm)
1491  else
1492 ! write(*,'(a,4i8)') 'different dimensions,nx,ny,nxc,nyc : ',nx,ny,nxc,nyc
1493  call intp_tile(tf_clm0, cxlats, cxlons, nyc, nxc, &
1494  tf_clm_ij,xlats_ij,xlons_ij,ny, nx)
1495  call intp_tile(tf_trd0, cxlats, cxlons, nyc, nxc, &
1496  tf_trd_ij,xlats_ij,xlons_ij,ny, nx)
1497 ! write(*,'(a,2f9.3)') 'tf_clm0, min, max : ',minval(tf_clm0),maxval(tf_clm0)
1498 
1499  tf_clm(:,:) = reshape(tf_clm_ij, (/nx,ny/) )
1500  tf_trd(:,:) = reshape(tf_trd_ij, (/nx,ny/) )
1501  endif
1502 
1503 end subroutine get_tf_clm
1504 
1519 subroutine get_tf_clm_ta(tf_clm_ta,tf_clm_trend,xlats,xlons,nlat,nlon,mon1,mon2,wei1,wei2)
1520  use read_write_data, only : read_tf_clim_grb
1521  implicit none
1522 
1523 ! input
1524  integer, intent(in) :: nlat,nlon,mon1,mon2
1525  real, intent(in) :: wei1,wei2
1526 ! output
1527  real, dimension(nlon,nlat), intent(out) :: tf_clm_ta,tf_clm_trend
1528  real, dimension(nlat), intent(out) :: xlats
1529  real, dimension(nlon), intent(out) :: xlons
1530 
1531 !input/output data file names
1532  character (len=6), parameter :: fin_tf_clm='sstclm'
1533 
1534 ! local declare
1535  real, dimension(nlon,nlat) :: tf_clm1,tf_clm2
1536 
1537 !
1538 ! read in rtg sst climatology without bitmap (surface mask) for mon1 and mon2
1539 !
1540  call read_tf_clim_grb(trim(fin_tf_clm),tf_clm1,xlats,xlons,nlat,nlon,mon1)
1541  call read_tf_clim_grb(trim(fin_tf_clm),tf_clm2,xlats,xlons,nlat,nlon,mon2)
1542 !
1543 ! tf_clim at the analysis time
1544 !
1545  tf_clm_ta(:,:) = wei1*tf_clm1(:,:)+wei2*tf_clm2(:,:)
1546 !
1547 ! tf_clim trend at the analysis time (month to month)
1548 !
1549  tf_clm_trend(:,:) = (tf_clm2(:,:)-tf_clm1(:,:))/120.0
1550 
1551  write(*,'(a,2f9.3)') 'tf_clm_ta, min, max : ',minval(tf_clm_ta),maxval(tf_clm_ta)
1552  write(*,'(a,2f9.3)') 'tf_clm_trend, min, max : ',minval(tf_clm_trend),maxval(tf_clm_trend)
1553  end subroutine get_tf_clm_ta
1554 
1567 subroutine get_sal_clm(xlats_ij,xlons_ij,ny,nx,iy,im,id,ih,sal_clm)
1568  use read_write_data, only : get_dim_nc
1569  implicit none
1570 
1571  real, dimension(nx*ny), intent(in) :: xlats_ij !
1572  real, dimension(nx*ny), intent(in) :: xlons_ij !
1573  real, dimension(nx,ny), intent(out) :: sal_clm !
1574  integer, intent(in) :: iy,im,id,ih,nx,ny
1575 ! local declare
1576  real, allocatable, dimension(:,:) :: sal_clm0 ! salinity climatology at the valid time
1577  real, allocatable, dimension(:) :: cxlats ! latitudes of sst climatology
1578  real, allocatable, dimension(:) :: cxlons ! longitudes of sst climatology
1579 
1580  real, dimension(nx*ny) :: sal_clm_ij ! salinity climatology at target grids (nx*ny)
1581  real :: wei1,wei2
1582  integer :: nxc,nyc,mon1,mon2,i,j
1583  character (len=6), parameter :: fin_sal_clm='salclm' ! salinity climatology file name
1584 !
1585 ! get which two months used and their weights from atime
1586 !
1587  call get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1588 !
1589 ! get the dimensions of the sst climatology & allocate the related arrays
1590 !
1591  call get_dim_nc(fin_sal_clm,nyc,nxc)
1592  allocate( sal_clm0(nxc,nyc),cxlats(nyc),cxlons(nxc) )
1593 !
1594 ! get sal_clm at the analysis time from monthly climatology & cxlats, cxlons
1595 !
1596  call get_sal_clm_ta(sal_clm0,cxlats,cxlons,nyc,nxc,mon1,mon2,wei1,wei2)
1597 !
1598 ! get sal_clm (nx by ny lat/lon) valid at atime
1599 !
1600  if ( nx == nxc .and. ny == nyc ) then
1601  sal_clm(:,:) = sal_clm0(:,:)
1602 ! write(*,'(a,2f9.3)') 'same dimensions, sal_clm, min : ',minval(sal_clm),maxval(sal_clm)
1603  else
1604 ! write(*,'(a,4i8)') 'different dimensions,nx,ny,nxc,nyc : ',nx,ny,nxc,nyc
1605  call intp_tile(sal_clm0, cxlats, cxlons, nyc, nxc, &
1606  sal_clm_ij,xlats_ij,xlons_ij,ny, nx)
1607 ! write(*,'(a,2f9.3)') 'sal_clm0, min, max : ',minval(sal_clm0),maxval(sal_clm0)
1608 ! write(*,'(a,2f9.3)') 'done with intp_tile for sal_clm, min, max : ',minval(sal_clm_ij),maxval(sal_clm_ij)
1609 
1610  sal_clm(:,:) = reshape(sal_clm_ij, (/nx,ny/) )
1611  endif
1612 
1613 end subroutine get_sal_clm
1614 
1627 subroutine get_sal_clm_ta(sal_clm_ta,xlats,xlons,nlat,nlon,mon1,mon2,wei1,wei2)
1628 
1629  use read_write_data, only : read_salclm_gfs_nc
1630  implicit none
1631 
1632 ! input
1633  integer, intent(in) :: nlat,nlon,mon1,mon2
1634  real, intent(in) :: wei1,wei2
1635 ! output
1636  real, dimension(nlon,nlat), intent(out) :: sal_clm_ta
1637  real, dimension(nlat), intent(out) :: xlats
1638  real, dimension(nlon), intent(out) :: xlons
1639 
1640 !input/output data file names
1641  character (len=6), parameter :: fin_sal_clm='salclm'
1642 
1643 ! local declare
1644  real, dimension(nlon,nlat) :: sal_clm1,sal_clm2
1645 
1646 !
1647 ! read in rtg sst climatology without bitmap (surface mask) for mon1 and mon2
1648 !
1649  call read_salclm_gfs_nc(trim(fin_sal_clm),sal_clm1,xlats,xlons,nlat,nlon,mon1)
1650  call read_salclm_gfs_nc(trim(fin_sal_clm),sal_clm2,xlats,xlons,nlat,nlon,mon2)
1651 !
1652 ! sal_clim at the analysis time
1653 !
1654  sal_clm_ta(:,:) = wei1*sal_clm1(:,:)+wei2*sal_clm2(:,:)
1655  write(*,'(a,2f9.3)') 'sal_clm_ta, min, max : ',minval(sal_clm_ta),maxval(sal_clm_ta)
1656  end subroutine get_sal_clm_ta
1657 
1672 subroutine intp_tile(tf_lalo,dlats_lalo,dlons_lalo,jdim_lalo,idim_lalo, &
1673  tf_tile,xlats_tile,xlons_tile,jdim_tile,idim_tile)
1674 
1675  implicit none
1676 
1677 ! input/output
1678  real, dimension(idim_lalo,jdim_lalo), intent(in) :: tf_lalo
1679  real, dimension(jdim_lalo), intent(in) :: dlats_lalo
1680  real, dimension(idim_lalo), intent(in) :: dlons_lalo
1681  real, dimension(jdim_tile*idim_tile), intent(in) :: xlats_tile
1682  real, dimension(jdim_tile*idim_tile), intent(in) :: xlons_tile
1683  integer, intent(in) :: jdim_lalo,idim_lalo,jdim_tile,idim_tile
1684  real, dimension(jdim_tile*idim_tile), intent(out) :: tf_tile
1685 
1686 ! local
1687  real, parameter :: deg2rad=3.1415926/180.0
1688  real, dimension(jdim_lalo) :: xlats_lalo
1689  real, dimension(idim_lalo) :: xlons_lalo
1690  real :: tf,wsum,res_km
1691  integer :: itile,jtile
1692  integer :: ii,jj,ij,iii,jjj
1693  integer :: ilalo,jlalo,ilalop1,jlalop1
1694  integer :: istart,iend,jstart,jend,krad
1695 
1696  integer, allocatable, dimension(:,:) :: id1,id2,jdc
1697  real, allocatable, dimension(:,:,:) :: agrid,s2c
1698 
1699  print*
1700  print*,'interpolate from lat/lon grids to any one grid with known lat/lon'
1701 
1702  xlats_lalo = dlats_lalo*deg2rad
1703  xlons_lalo = dlons_lalo*deg2rad
1704 
1705  allocate(agrid(idim_tile,jdim_tile,2))
1706  agrid(:,:,1) = reshape(xlons_tile, (/idim_tile,jdim_tile/) )
1707  agrid(:,:,2) = reshape(xlats_tile, (/idim_tile,jdim_tile/) )
1708  agrid = agrid*deg2rad
1709 
1710  allocate(id1(idim_tile,jdim_tile))
1711  allocate(id2(idim_tile,jdim_tile))
1712  allocate(jdc(idim_tile,jdim_tile))
1713  allocate(s2c(idim_tile,jdim_tile,4))
1714 
1715 !----------------------------------------------------------------------
1716 ! compute bilinear weights for each model point from the nearest
1717 ! four lalo points. does not account for mask. that
1718 ! happens later.
1719 !----------------------------------------------------------------------
1720 
1721  call remap_coef( 1, idim_tile, 1, jdim_tile, idim_lalo, jdim_lalo, &
1722  xlons_lalo, xlats_lalo, id1, id2, jdc, s2c, agrid )
1723 
1724  do ij = 1, jdim_tile*idim_tile
1725 
1726  jtile = (ij-1)/idim_tile + 1
1727  itile = mod(ij,idim_tile)
1728  if (itile==0) itile = idim_tile
1729 
1730  ilalo = id1(itile,jtile)
1731  ilalop1 = id2(itile,jtile)
1732  jlalo = jdc(itile,jtile)
1733  jlalop1 = jdc(itile,jtile) + 1
1734 
1735  wsum = s2c(itile,jtile,1) + s2c(itile,jtile,2) + &
1736  s2c(itile,jtile,3) + s2c(itile,jtile,4)
1737 
1738  tf_tile(ij) = ( s2c(itile,jtile,1)*tf_lalo(ilalo,jlalo) + &
1739  s2c(itile,jtile,2)*tf_lalo(ilalop1,jlalo) + &
1740  s2c(itile,jtile,3)*tf_lalo(ilalop1,jlalop1) + &
1741  s2c(itile,jtile,4)*tf_lalo(ilalo,jlalop1) )/wsum
1742  enddo
1743 
1744  deallocate(id1, id2, jdc, s2c)
1745 
1746 end subroutine intp_tile
1747 
1760 subroutine get_tim_wei(iy,im,id,ih,mon1,mon2,wei1,wei2)
1761  implicit none
1762 
1763 ! input
1764  integer, intent(in) :: iy,im,id,ih
1765 ! output
1766  integer, intent(out) :: mon1,mon2
1767  real, intent(out) :: wei1,wei2
1768 
1769 ! local declare
1770  real :: rjday
1771  integer :: mon,monend,monm,monp,jdow,jdoy,jday
1772  integer :: jda(8)
1773 !
1774 !dayhf : julian day of the middle of each month
1775 !
1776  real, dimension(13) :: dayhf
1777  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/
1778 
1779 ! 15, 44, 73.5, 104, 134.5, 165, 195.5, 226.5, 257, 287.5, 318.5, 349 ! from
1780 ! woa05
1781 
1782  jda=0
1783  jda(1)=iy
1784  jda(2)=im
1785  jda(3)=id
1786  jda(5)=ih
1787  jdow = 0
1788  jdoy = 0
1789  jday = 0
1790  call w3doxdat(jda,jdow,jdoy,jday)
1791  rjday=jdoy+jda(5)/24.
1792  if(rjday.lt.dayhf(1)) rjday=rjday+365.
1793 
1794  monend = 12
1795  do mon = 1, monend
1796  monm = mon
1797  monp = mon + 1
1798  if( rjday >= dayhf(monm) .and. rjday < dayhf(monp) ) then
1799  mon1 = monm
1800  mon2 = monp
1801  go to 10
1802  endif
1803  enddo
1804 
1805  print *,'FATAL ERROR in get_tim_wei, wrong rjday',rjday
1806  call abort
1807  10 continue
1808 
1809  wei1 = (dayhf(mon2)-rjday)/(dayhf(mon2)-dayhf(mon1))
1810  wei2 = (rjday-dayhf(mon1))/(dayhf(mon2)-dayhf(mon1))
1811 
1812  if( mon2 == 13 ) mon2=1
1813 
1814  write(*,'(a,2i4,3f9.3)') 'mon1,mon2,rjday,wei1,wei2=',mon1,mon2,rjday,wei1,wei2
1815 
1816  end subroutine get_tim_wei
1817 
1827  real function tfreez(salinity)
1828 
1829  implicit none
1830 
1831  REAL salinity,sal
1832  REAL a1, a2, a3
1833  parameter(a1 = -0.0575)
1834  parameter(a2 = 1.710523e-3)
1835  parameter(a3 = -2.154996e-4)
1836 
1837  IF (salinity .LT. 0.) THEN
1838  sal = 0.
1839  ELSE
1840  sal = salinity
1841  ENDIF
1842  tfreez = sal*(a1+a2*sqrt(sal)+a3*sal)
1843 
1844  return
1845  end
1846 
num_parthds
integer function num_parthds()
Return the number of omp threads.
Definition: num_parthds.f90:9
adjust_nsst
subroutine adjust_nsst(RLA, RLO, SLMSK_TILE, SLMSK_FG_TILE, SKINT_TILE, SICET_TILE, sice_tile, SOILT_TILE, NSST, LENSFC, LSOIL, IDIM, JDIM, ZSEA1, ZSEA2, MON, DAY, DELTSFC, tf_clm_tile, tf_trd_tile, sal_clm_tile)
Read in gsi file with the updated reference temperature increments (on the gaussian grid)...
Definition: cycle.f90:575
read_write_data::read_gsi_data
subroutine, public read_gsi_data(GSI_FILE)
Read file from the GSI containing the foundation temperature increments and mask. ...
Definition: read_write_data.f90:1075
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:1450
read_write_data::read_lat_lon_orog
subroutine, public read_lat_lon_orog(RLA, RLO, OROG, OROG_UF, TILE_NUM, IDIM, JDIM, IJDIM)
Read latitude and longitude for the cubed-sphere tile from the &#39;grid&#39; file.
Definition: read_write_data.f90:911
read_write_data::read_data
subroutine, public read_data(TSFFCS, SMCFCS, SNOFCS, STCFCS, TG3FCS, ZORFCS, CVFCS, CVBFCS, CVTFCS, ALBFCS, SLIFCS, VEGFCS, CNPFCS, F10M, VETFCS, SOTFCS, ALFFCS, USTAR, FMM, FHH, SIHFCS, SICFCS, SITFCS, TPRCP, SRFLAG, SWDFCS, VMNFCS, VMXFCS, SLCFCS, SLPFCS, ABSFCS, T2M, Q2M, SLMASK, ZSOIL, LSOIL, LENSFC, DO_NSST, NSST)
Read the first guess surface records and nsst records (if selected) for a single cubed-sphere tile...
Definition: read_write_data.f90:1180
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:1519
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:1136
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:1627
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:1313
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:1897
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:1821
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:1672
remap_coef
subroutine 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: cycle.f90:1214
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:983
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:1762
read_write_data::write_data
subroutine, public write_data(slifcs, tsffcs, snofcs, tg3fcs, zorfcs, albfcs, alffcs, vegfcs, cnpfcs, f10m, t2m, q2m, vetfcs, sotfcs, ustar, fmm, fhh, sicfcs, sihfcs, sitfcs, tprcp, srflag, swdfcs, vmnfcs, vmxfcs, slpfcs, absfcs, slcfcs, smcfcs, stcfcs, idim, jdim, lensfc, lsoil, do_nsst, nsst)
Write out all surface records - and nsst records if selected - on a single cubed-sphere tile to a mod...
Definition: read_write_data.f90:108
sfc_drv
program sfc_drv
Stand alone surface/NSST cycle driver for the cubed-sphere grid.
Definition: cycle.f90:89
sfccycle
subroutine sfccycle(lugb, len, lsoil, sig1t, deltsfc, iy, im, id, ih, fh, rla, rlo, slmask, orog, orog_uf, use_ufo, nst_anl, sihfcs, sicfcs, sitfcs, swdfcs, slcfcs, vmnfcs, vmxfcs, slpfcs, absfcs, tsffcs, snofcs, zorfcs, albfcs, tg3fcs, cnpfcs, smcfcs, stcfcs, slifcs, aisfcs, vegfcs, vetfcs, sotfcs, alffcs, cvfcs, cvbfcs, cvtfcs, me, nlunit, sz_nml, input_nml_file, ialb, isot, ivegsrc, tile_num_ch, i_index, j_index)
Surface cycling driver routine.
Definition: sfcsub.F:175
sfcdrv
subroutine sfcdrv(LUGB, IDIM, JDIM, LENSFC, LSOIL, DELTSFC, IY, IM, ID, IH, FH, IALB, USE_UFO, DO_NSST, ADJT_NST_ONLY, ZSEA1, ZSEA2, ISOT, IVEGSRC, MYRANK)
Driver routine for updating the surface fields.
Definition: cycle.f90:283
tfreez
real function tfreez(salinity)
Compute the freezing point of water as a function of salinity.
Definition: cycle.f90:1827
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:1760
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:1567
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:1404
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:1620