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