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