New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 13540 for NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC – NEMO

Ignore:
Timestamp:
2020-09-29T12:41:06+02:00 (4 years ago)
Author:
andmirek
Message:

Ticket #2386: update to latest trunk

Location:
NEMO/branches/2020/r12377_ticket2386
Files:
27 edited

Legend:

Unmodified
Added
Removed

  • NEMO/branches/2020/r12377_ticket2386

    • Property svn:externals

      •  

        old new  
        33^/utils/build/mk@HEAD         mk 
        44^/utils/tools@HEAD            tools 
        5 ^/vendors/AGRIF/dev@HEAD      ext/AGRIF 
         5^/vendors/AGRIF/dev_r12970_AGRIF_CMEMS      ext/AGRIF 
        66^/vendors/FCM@HEAD            ext/FCM 
        77^/vendors/IOIPSL@HEAD         ext/IOIPSL 
        88 
        99# SETTE 
        10 ^/utils/CI/sette@HEAD         sette 
         10^/utils/CI/sette@13507        sette 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/cpl_oasis3.F90

    r12377 r13540  
    6969   INTEGER, PUBLIC, PARAMETER ::   nmaxcat=5    ! Maximum number of coupling fields 
    7070   INTEGER, PUBLIC, PARAMETER ::   nmaxcpl=5    ! Maximum number of coupling fields 
    71    LOGICAL, PARAMETER         ::   ltmp_wapatch = .TRUE.   ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define   
    72    INTEGER                    ::   nldi_save, nlei_save 
    73    INTEGER                    ::   nldj_save, nlej_save 
    7471    
    7572   TYPE, PUBLIC ::   FLD_CPL               !: Type for coupling field information 
    … …  
    148145      !!-------------------------------------------------------------------- 
    149146 
    150       ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 
    151       IF( ltmp_wapatch ) THEN 
    152          nldi_save = nldi   ;   nlei_save = nlei 
    153          nldj_save = nldj   ;   nlej_save = nlej 
    154          IF( nimpp           ==      1 ) nldi = 1 
    155          IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi 
    156          IF( njmpp           ==      1 ) nldj = 1 
    157          IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj 
    158       ENDIF  
    159147      IF(lwp) WRITE(numout,*) 
    160148      IF(lwp) WRITE(numout,*) 'cpl_define : initialization in coupled ocean/atmosphere case' 
    … …  
    177165      ENDIF 
    178166      ! 
    179       ! ... Define the shape for the area that excludes the halo 
    180       !     For serial configuration (key_mpp_mpi not being active) 
    181       !     nl* is set to the global values 1 and jp*glo. 
     167      ! ... Define the shape for the area that excludes the halo as we don't want them to be "seen" by oasis 
    182168      ! 
    183169      ishape(1) = 1 
    184       ishape(2) = nlei-nldi+1 
     170      ishape(2) = Ni_0 
    185171      ishape(3) = 1 
    186       ishape(4) = nlej-nldj+1 
     172      ishape(4) = Nj_0 
    187173      ! 
    188174      ! ... Allocate memory for data exchange 
    189175      ! 
    190       ALLOCATE(exfld(nlei-nldi+1, nlej-nldj+1), stat = nerror) 
     176      ALLOCATE(exfld(Ni_0, Nj_0), stat = nerror)        ! allocate only inner domain (without halos) 
    191177      IF( nerror > 0 ) THEN 
    192178         CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in allocating exfld')   ;   RETURN 
    … …  
    194180      ! 
    195181      ! ----------------------------------------------------------------- 
    196       ! ... Define the partition  
     182      ! ... Define the partition, excluding halos as we don't want them to be "seen" by oasis     
    197183      ! ----------------------------------------------------------------- 
    198184       
    199       paral(1) = 2                                              ! box partitioning 
    200       paral(2) = jpiglo * (nldj-1+njmpp-1) + (nldi-1+nimpp-1)   ! NEMO lower left corner global offset     
    201       paral(3) = nlei-nldi+1                                    ! local extent in i  
    202       paral(4) = nlej-nldj+1                                    ! local extent in j 
    203       paral(5) = jpiglo                                         ! global extent in x 
     185      paral(1) = 2                                      ! box partitioning 
     186      paral(2) = Ni0glo * mjg0(nn_hls) + mig0(nn_hls)   ! NEMO lower left corner global offset, without halos  
     187      paral(3) = Ni_0                                   ! local extent in i, excluding halos 
     188      paral(4) = Nj_0                                   ! local extent in j, excluding halos 
     189      paral(5) = Ni0glo                                 ! global extent in x, excluding halos 
    204190       
    205191      IF( sn_cfctl%l_oasout ) THEN 
    206192         WRITE(numout,*) ' multiexchg: paral (1:5)', paral 
    207          WRITE(numout,*) ' multiexchg: jpi, jpj =', jpi, jpj 
    208          WRITE(numout,*) ' multiexchg: nldi, nlei, nimpp =', nldi, nlei, nimpp 
    209          WRITE(numout,*) ' multiexchg: nldj, nlej, njmpp =', nldj, nlej, njmpp 
     193         WRITE(numout,*) ' multiexchg: Ni_0, Nj_0 =', Ni_0, Nj_0 
     194         WRITE(numout,*) ' multiexchg: Nis0, Nie0, nimpp =', Nis0, Nie0, nimpp 
     195         WRITE(numout,*) ' multiexchg: Njs0, Nje0, njmpp =', Njs0, Nje0, njmpp 
    210196      ENDIF 
    211197    
    212       CALL oasis_def_partition ( id_part, paral, nerror, jpiglo*jpjglo ) 
     198      CALL oasis_def_partition ( id_part, paral, nerror, Ni0glo*Nj0glo )   ! global number of points, excluding halos 
    213199      ! 
    214200      ! ... Announce send variables.  
    … …  
    316302#endif 
    317303      ! 
    318       IF( ltmp_wapatch ) THEN 
    319          nldi = nldi_save   ;   nlei = nlei_save 
    320          nldj = nldj_save   ;   nlej = nlej_save 
    321       ENDIF 
    322304   END SUBROUTINE cpl_define 
    323305    
    … …  
    337319      INTEGER                                   ::   jc,jm     ! local loop index 
    338320      !!-------------------------------------------------------------------- 
    339       ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 
    340       IF( ltmp_wapatch ) THEN 
    341          nldi_save = nldi   ;   nlei_save = nlei 
    342          nldj_save = nldj   ;   nlej_save = nlej 
    343          IF( nimpp           ==      1 ) nldi = 1 
    344          IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi 
    345          IF( njmpp           ==      1 ) nldj = 1 
    346          IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj 
    347       ENDIF 
    348321      ! 
    349322      ! snd data to OASIS3 
    … …  
    352325         DO jm = 1, ssnd(kid)%ncplmodel 
    353326         
    354             IF( ssnd(kid)%nid(jc,jm) /= -1 ) THEN 
    355                CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(nldi:nlei, nldj:nlej,jc), kinfo ) 
     327            IF( ssnd(kid)%nid(jc,jm) /= -1 ) THEN   ! exclude halos from data sent to oasis 
     328               CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(Nis0:Nie0, Njs0:Nje0,jc), kinfo ) 
    356329                
    357330               IF ( sn_cfctl%l_oasout ) THEN         
    … …  
    363336                     WRITE(numout,*) 'oasis_put:  kstep ', kstep 
    364337                     WRITE(numout,*) 'oasis_put:   info ', kinfo 
    365                      WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(nldi:nlei,nldj:nlej,jc)) 
    366                      WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(nldi:nlei,nldj:nlej,jc)) 
    367                      WRITE(numout,*) '     -     Sum value is ', SUM(pdata(nldi:nlei,nldj:nlej,jc)) 
     338                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc)) 
     339                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc)) 
     340                     WRITE(numout,*) '     -     Sum value is ',    SUM(pdata(Nis0:Nie0,Njs0:Nje0,jc)) 
    368341                     WRITE(numout,*) '****************' 
    369342                  ENDIF 
    … …  
    374347         ENDDO 
    375348      ENDDO 
    376       IF( ltmp_wapatch ) THEN 
    377          nldi = nldi_save   ;   nlei = nlei_save 
    378          nldj = nldj_save   ;   nlej = nlej_save 
    379       ENDIF 
    380349      ! 
    381350    END SUBROUTINE cpl_snd 
    … …  
    396365      !! 
    397366      INTEGER                                   ::   jc,jm     ! local loop index 
    398       LOGICAL                                   ::   llaction, llfisrt 
     367      LOGICAL                                   ::   llaction, ll_1st 
    399368      !!-------------------------------------------------------------------- 
    400       ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 
    401       IF( ltmp_wapatch ) THEN 
    402          nldi_save = nldi   ;   nlei_save = nlei 
    403          nldj_save = nldj   ;   nlej_save = nlej 
    404       ENDIF 
    405369      ! 
    406370      ! receive local data from OASIS3 on every process 
    … …  
    409373      ! 
    410374      DO jc = 1, srcv(kid)%nct 
    411          IF( ltmp_wapatch ) THEN 
    412             IF( nimpp           ==      1 ) nldi = 1 
    413             IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi 
    414             IF( njmpp           ==      1 ) nldj = 1 
    415             IF( njmpp + jpj - 1 == jpjglo ) nlej = jpj 
    416          ENDIF 
    417          llfisrt = .TRUE. 
     375         ll_1st = .TRUE. 
    418376 
    419377         DO jm = 1, srcv(kid)%ncplmodel 
    … …  
    426384                  &        kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut 
    427385                
    428                IF ( sn_cfctl%l_oasout )   WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm) 
     386               IF ( sn_cfctl%l_oasout )   & 
     387                  &  WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm) 
    429388                
    430                IF( llaction ) THEN 
     389               IF( llaction ) THEN   ! data received from oasis do not include halos 
    431390                   
    432391                  kinfo = OASIS_Rcv 
    433                   IF( llfisrt ) THEN  
    434                      pdata(nldi:nlei,nldj:nlej,jc) =                                 exfld(:,:) * pmask(nldi:nlei,nldj:nlej,jm) 
    435                      llfisrt = .FALSE. 
     392                  IF( ll_1st ) THEN  
     393                     pdata(Nis0:Nie0,Njs0:Nje0,jc) =   exfld(:,:) * pmask(Nis0:Nie0,Njs0:Nje0,jm) 
     394                     ll_1st = .FALSE. 
    436395                  ELSE 
    437                      pdata(nldi:nlei,nldj:nlej,jc) = pdata(nldi:nlei,nldj:nlej,jc) + exfld(:,:) * pmask(nldi:nlei,nldj:nlej,jm) 
     396                     pdata(Nis0:Nie0,Njs0:Nje0,jc) = pdata(Nis0:Nie0,Njs0:Nje0,jc)   & 
     397                        &                                + exfld(:,:) * pmask(Nis0:Nie0,Njs0:Nje0,jm) 
    438398                  ENDIF 
    439399                   
    … …  
    444404                     WRITE(numout,*) 'oasis_get:   kstep', kstep 
    445405                     WRITE(numout,*) 'oasis_get:   info ', kinfo 
    446                      WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(:,:,jc)) 
    447                      WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(:,:,jc)) 
    448                      WRITE(numout,*) '     -     Sum value is ', SUM(pdata(:,:,jc)) 
     406                     WRITE(numout,*) '     - Minimum value is ', MINVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc)) 
     407                     WRITE(numout,*) '     - Maximum value is ', MAXVAL(pdata(Nis0:Nie0,Njs0:Nje0,jc)) 
     408                     WRITE(numout,*) '     -     Sum value is ',    SUM(pdata(Nis0:Nie0,Njs0:Nje0,jc)) 
    449409                     WRITE(numout,*) '****************' 
    450410                  ENDIF 
    … …  
    456416         ENDDO 
    457417 
    458          IF( ltmp_wapatch ) THEN 
    459             nldi = nldi_save   ;   nlei = nlei_save 
    460             nldj = nldj_save   ;   nlej = nlej_save 
    461          ENDIF 
    462          !--- Fill the overlap areas and extra hallows (mpp) 
    463          !--- check periodicity conditions (all cases) 
    464          IF( .not. llfisrt ) THEN 
     418         !--- we must call lbc_lnk to fill the halos that where not received. 
     419         IF( .NOT. ll_1st ) THEN 
    465420            CALL lbc_lnk( 'cpl_oasis3', pdata(:,:,jc), srcv(kid)%clgrid, srcv(kid)%nsgn )    
    466421         ENDIF 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/cyclone.F90

    r12377 r13540  
    147147            zb = 2. 
    148148 
    149             DO_2D_11_11 
     149            DO_2D( 1, 1, 1, 1 ) 
    150150 
    151151               ! calc distance between TC center and any point following great circle 
    … …  
    208208            ENDIF            
    209209         
    210             DO_2D_11_11 
     210            DO_2D( 1, 1, 1, 1 ) 
    211211                                
    212212               zzrglam = rad * glamt(ji,jj) - zrlon 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/fldread.F90

    r12511 r13540  
    5353      LOGICAL              ::   ln_tint     ! time interpolation or not (T/F) 
    5454      LOGICAL              ::   ln_clim     ! climatology or not (T/F) 
    55       CHARACTER(len = 8)   ::   cltype      ! type of data file 'daily', 'monthly' or yearly' 
     55      CHARACTER(len = 8)   ::   clftyp      ! type of data file 'daily', 'monthly' or yearly' 
    5656      CHARACTER(len = 256) ::   wname       ! generic name of a NetCDF weights file to be used, blank if not 
    5757      CHARACTER(len = 34)  ::   vcomp       ! symbolic component name if a vector that needs rotation 
    … …  
    6969      LOGICAL                         ::   ln_tint      ! time interpolation or not (T/F) 
    7070      LOGICAL                         ::   ln_clim      ! climatology or not (T/F) 
    71       CHARACTER(len = 8)              ::   cltype       ! type of data file 'daily', 'monthly' or yearly' 
     71      CHARACTER(len = 8)              ::   clftyp       ! type of data file 'daily', 'monthly' or yearly' 
     72      CHARACTER(len = 1)              ::   cltype       ! nature of grid-points: T, U, V... 
     73      REAL(wp)                        ::   zsgn         ! -1. the sign change across the north fold, =  1. otherwise 
    7274      INTEGER                         ::   num          ! iom id of the jpfld files to be read 
    73       INTEGER , DIMENSION(2)          ::   nrec_b       ! before record (1: index, 2: second since Jan. 1st 00h of nit000 year) 
    74       INTEGER , DIMENSION(2)          ::   nrec_a       ! after  record (1: index, 2: second since Jan. 1st 00h of nit000 year) 
    75       INTEGER , ALLOCATABLE, DIMENSION(:      ) ::   nrecsec   !  
    76       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:  ) ::   fnow   ! input fields interpolated to now time step 
    77       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   fdta   ! 2 consecutive record of input fields 
     75      INTEGER , DIMENSION(2,2)        ::   nrec         ! before/after record (1: index, 2: second since Jan. 1st 00h of yr nit000) 
     76      INTEGER                         ::   nbb          ! index of before values 
     77      INTEGER                         ::   naa          ! index of after  values 
     78      INTEGER , ALLOCATABLE, DIMENSION(:) ::   nrecsec   !  
     79      REAL(wp), POINTER, DIMENSION(:,:,:  ) ::   fnow   ! input fields interpolated to now time step 
     80      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   fdta   ! 2 consecutive record of input fields 
    7881      CHARACTER(len = 256)            ::   wgtname      ! current name of the NetCDF weight file acting as a key 
    7982      !                                                 ! into the WGTLIST structure 
    … …  
    127130   !! * Substitutions 
    128131#  include "do_loop_substitute.h90" 
     132#  include "domzgr_substitute.h90" 
    129133   !!---------------------------------------------------------------------- 
    130134   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    156160      INTEGER  ::   jf           ! dummy indices 
    157161      INTEGER  ::   isecsbc      ! number of seconds between Jan. 1st 00h of nit000 year and the middle of sbc time step 
     162      INTEGER  ::   ibb, iaa     ! shorter name for sd(jf)%nbb and sd(jf)%naa 
    158163      LOGICAL  ::   ll_firstcall ! true if this is the first call to fld_read for this set of fields 
    159164      REAL(wp) ::   zt_offset    ! local time offset variable 
    … …  
    203208            IF( TRIM(sd(jf)%clrootname) == 'NOT USED' )   CYCLE 
    204209            ! 
     210            ibb = sd(jf)%nbb   ;   iaa = sd(jf)%naa 
     211            ! 
    205212            IF( sd(jf)%ln_tint ) THEN              ! temporal interpolation 
    206213               IF(lwp .AND. kt - nit000 <= 100 ) THEN  
    … …  
    208215                     &    "', records b/a: ', i6.4, '/', i6.4, ' (days ', f9.4,'/', f9.4, ')')" 
    209216                  WRITE(numout, clfmt)  TRIM( sd(jf)%clvar ), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,   &             
    210                      & sd(jf)%nrec_b(1), sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday 
     217                     & sd(jf)%nrec(1,ibb), sd(jf)%nrec(1,iaa), REAL(sd(jf)%nrec(2,ibb),wp)/rday, REAL(sd(jf)%nrec(2,iaa),wp)/rday 
    211218                  WRITE(numout, *) '      zt_offset is : ',zt_offset 
    212219               ENDIF 
    213220               ! temporal interpolation weights 
    214                ztinta =  REAL( isecsbc - sd(jf)%nrec_b(2), wp ) / REAL( sd(jf)%nrec_a(2) - sd(jf)%nrec_b(2), wp ) 
     221               ztinta =  REAL( isecsbc - sd(jf)%nrec(2,ibb), wp ) / REAL( sd(jf)%nrec(2,iaa) - sd(jf)%nrec(2,ibb), wp ) 
    215222               ztintb =  1. - ztinta 
    216                sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,1) + ztinta * sd(jf)%fdta(:,:,:,2) 
     223               sd(jf)%fnow(:,:,:) = ztintb * sd(jf)%fdta(:,:,:,ibb) + ztinta * sd(jf)%fdta(:,:,:,iaa) 
    217224            ELSE   ! nothing to do... 
    218225               IF(lwp .AND. kt - nit000 <= 100 ) THEN 
    … …  
    220227                     &    "', record: ', i6.4, ' (days ', f9.4, ' <-> ', f9.4, ')')" 
    221228                  WRITE(numout, clfmt) TRIM(sd(jf)%clvar), kt, REAL(isecsbc,wp)/rday, nyear, nmonth, nday,    & 
    222                      &                 sd(jf)%nrec_a(1), REAL(sd(jf)%nrec_b(2),wp)/rday, REAL(sd(jf)%nrec_a(2),wp)/rday 
     229                     &                 sd(jf)%nrec(1,iaa), REAL(sd(jf)%nrec(2,ibb),wp)/rday, REAL(sd(jf)%nrec(2,iaa),wp)/rday 
    223230               ENDIF 
    224231            ENDIF 
    … …  
    250257      ! 
    251258      CALL fld_clopn( sdjf ) 
    252       sdjf%nrec_a(:) = (/ 1, nflag /)  ! default definition to force flp_update to read the file. 
     259      sdjf%nrec(:,sdjf%naa) = (/ 1, nflag /)  ! default definition to force flp_update to read the file. 
    253260      ! 
    254261   END SUBROUTINE fld_init 
    … …  
    261268      !! ** Purpose : Compute 
    262269      !!              if sdjf%ln_tint = .TRUE. 
    263       !!                  nrec_a: record number and its time (nrec_b is obtained from nrec_a when swapping) 
     270      !!                  nrec(:,iaa): record number and its time (nrec(:,ibb) is obtained from nrec(:,iaa) when swapping) 
    264271      !!              if sdjf%ln_tint = .FALSE. 
    265       !!                  nrec_a(1): record number 
    266       !!                  nrec_b(2) and nrec_a(2): time of the beginning and end of the record 
     272      !!                  nrec(1,iaa): record number 
     273      !!                  nrec(2,ibb) and nrec(2,iaa): time of the beginning and end of the record 
    267274      !!---------------------------------------------------------------------- 
    268275      INTEGER  ,           INTENT(in   ) ::   ksecsbc   !  
    … …  
    270277      INTEGER  , OPTIONAL, INTENT(in   ) ::   Kmm    ! ocean time level index 
    271278      ! 
    272       INTEGER  ::   ja     ! end of this record (in seconds) 
    273       !!---------------------------------------------------------------------- 
    274       ! 
    275       IF( ksecsbc > sdjf%nrec_a(2) ) THEN     ! --> we need to update after data 
     279      INTEGER  ::   ja           ! end of this record (in seconds) 
     280      INTEGER  ::   ibb, iaa     ! shorter name for sdjf%nbb and sdjf%naa 
     281      !!---------------------------------------------------------------------- 
     282      ibb = sdjf%nbb   ;   iaa = sdjf%naa 
     283      ! 
     284      IF( ksecsbc > sdjf%nrec(2,iaa) ) THEN     ! --> we need to update after data 
    276285         
    277          ! find where is the new after record... (it is not necessary sdjf%nrec_a(1)+1 ) 
    278          ja = sdjf%nrec_a(1) 
     286         ! find where is the new after record... (it is not necessary sdjf%nrec(1,iaa)+1 ) 
     287         ja = sdjf%nrec(1,iaa) 
    279288         DO WHILE ( ksecsbc >= sdjf%nrecsec(ja) .AND. ja < sdjf%nreclast )   ! Warning: make sure ja <= sdjf%nreclast in this test 
    280289            ja = ja + 1 
    … …  
    283292 
    284293         ! if ln_tint and if the new after is not ja+1, we need also to update after data before the swap 
    285          ! so, after the swap, sdjf%nrec_b(2) will still be the closest value located just before ksecsbc 
    286          IF( sdjf%ln_tint .AND. ( ja > sdjf%nrec_a(1) + 1 .OR. sdjf%nrec_a(2) == nflag ) ) THEN 
    287             sdjf%nrec_a(:) = (/ ja-1, sdjf%nrecsec(ja-1) /)   ! update nrec_a with before information 
    288             CALL fld_get( sdjf, Kmm )                         ! read after data that will be used as before data 
     294         ! so, after the swap, sdjf%nrec(2,ibb) will still be the closest value located just before ksecsbc 
     295         IF( sdjf%ln_tint .AND. ( ja > sdjf%nrec(1,iaa) + 1 .OR. sdjf%nrec(2,iaa) == nflag ) ) THEN 
     296            sdjf%nrec(:,iaa) = (/ ja-1, sdjf%nrecsec(ja-1) /)   ! update nrec(:,iaa) with before information 
     297            CALL fld_get( sdjf, Kmm )                           ! read after data that will be used as before data 
    289298         ENDIF 
    290299             
    … …  
    309318            ! if ln_tint and if after is not the first record, we must (potentially again) update after data before the swap 
    310319            IF( sdjf%ln_tint .AND. ja > 1 ) THEN 
    311                IF( sdjf%nrecsec(0) /= nflag ) THEN                  ! no trick used: after file is not the current file 
    312                   sdjf%nrec_a(:) = (/ ja-1, sdjf%nrecsec(ja-1) /)   ! update nrec_a with before information 
    313                   CALL fld_get( sdjf, Kmm )                         ! read after data that will be used as before data 
     320               IF( sdjf%nrecsec(0) /= nflag ) THEN                    ! no trick used: after file is not the current file 
     321                  sdjf%nrec(:,iaa) = (/ ja-1, sdjf%nrecsec(ja-1) /)   ! update nrec(:,iaa) with before information 
     322                  CALL fld_get( sdjf, Kmm )                           ! read after data that will be used as before data 
    314323               ENDIF 
    315324            ENDIF 
    … …  
    317326         ENDIF 
    318327 
    319          IF( sdjf%ln_tint ) THEN  
    320             ! Swap data 
    321             sdjf%nrec_b(:)     = sdjf%nrec_a(:)                     ! swap before record informations 
    322             sdjf%rotn(1)       = sdjf%rotn(2)                       ! swap before rotate informations 
    323             sdjf%fdta(:,:,:,1) = sdjf%fdta(:,:,:,2)                 ! swap before record field 
    324          ELSE 
    325             sdjf%nrec_b(:) = (/ ja-1, sdjf%nrecsec(ja-1) /)         ! only for print  
     328         IF( sdjf%ln_tint ) THEN                                ! Swap data 
     329            sdjf%nbb = sdjf%naa                                 !    swap indices 
     330            sdjf%naa = 3 - sdjf%naa                             !    = 2(1) if naa == 1(2) 
     331         ELSE                                                   ! No swap 
     332            sdjf%nrec(:,ibb) = (/ ja-1, sdjf%nrecsec(ja-1) /)   !    only for print  
    326333         ENDIF 
    327334             
    328335         ! read new after data 
    329          sdjf%nrec_a(:) = (/ ja, sdjf%nrecsec(ja) /)                ! update nrec_a as it is used by fld_get 
    330          CALL fld_get( sdjf, Kmm )                                  ! read after data (with nrec_a informations) 
     336         sdjf%nrec(:,sdjf%naa) = (/ ja, sdjf%nrecsec(ja) /)     ! update nrec(:,naa) as it is used by fld_get 
     337         CALL fld_get( sdjf, Kmm )                              ! read after data (with nrec(:,naa) informations) 
    331338         
    332339      ENDIF 
    … …  
    345352      ! 
    346353      INTEGER ::   ipk      ! number of vertical levels of sdjf%fdta ( 2D: ipk=1 ; 3D: ipk=jpk ) 
     354      INTEGER ::   iaa      ! shorter name for sdjf%naa 
    347355      INTEGER ::   iw       ! index into wgts array 
    348       INTEGER ::   ipdom    ! index of the domain 
    349356      INTEGER ::   idvar    ! variable ID 
    350357      INTEGER ::   idmspc   ! number of spatial dimensions 
    351358      LOGICAL ::   lmoor    ! C1D case: point data 
    352       !!--------------------------------------------------------------------- 
    353       ! 
    354       ipk = SIZE( sdjf%fnow, 3 ) 
    355       ! 
    356       IF( ASSOCIATED(sdjf%imap) ) THEN 
    357          IF( sdjf%ln_tint ) THEN   ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1),   & 
    358             &                                        sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint, Kmm ) 
    359          ELSE                      ;   CALL fld_map( sdjf%num, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1),   & 
    360             &                                        sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint, Kmm ) 
    361          ENDIF 
    362       ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN 
     359      REAL(wp), DIMENSION(:,:,:), POINTER ::   dta_alias   ! short cut 
     360      !!--------------------------------------------------------------------- 
     361      iaa = sdjf%naa 
     362      ! 
     363      IF( sdjf%ln_tint ) THEN   ;   dta_alias => sdjf%fdta(:,:,:,iaa) 
     364      ELSE                      ;   dta_alias => sdjf%fnow(:,:,:    ) 
     365      ENDIF 
     366      ipk = SIZE( dta_alias, 3 ) 
     367      ! 
     368      IF( ASSOCIATED(sdjf%imap) ) THEN              ! BDY case  
     369         CALL fld_map( sdjf%num, sdjf%clvar, dta_alias(:,:,:), sdjf%nrec(1,iaa),   & 
     370            &          sdjf%imap, sdjf%igrd, sdjf%ibdy, sdjf%ltotvel, sdjf%lzint, Kmm ) 
     371      ELSE IF( LEN(TRIM(sdjf%wgtname)) > 0 ) THEN   ! On-the-fly interpolation 
    363372         CALL wgt_list( sdjf, iw ) 
    364          IF( sdjf%ln_tint ) THEN   ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fdta(:,:,:,2),          &  
    365             &                                                                          sdjf%nrec_a(1), sdjf%lsmname ) 
    366          ELSE                      ;   CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, sdjf%fnow(:,:,:  ),          & 
    367             &                                                                          sdjf%nrec_a(1), sdjf%lsmname ) 
    368          ENDIF 
    369       ELSE 
    370          IF( SIZE(sdjf%fnow, 1) == jpi ) THEN   ;   ipdom = jpdom_data 
    371          ELSE                                   ;   ipdom = jpdom_unknown 
    372          ENDIF 
     373         CALL fld_interp( sdjf%num, sdjf%clvar, iw, ipk, dta_alias(:,:,:), sdjf%nrec(1,iaa), sdjf%lsmname ) 
     374         CALL lbc_lnk( 'fldread', dta_alias(:,:,:), sdjf%cltype, sdjf%zsgn, kfillmode = jpfillcopy ) 
     375      ELSE                                          ! default case 
    373376         ! C1D case: If product of spatial dimensions == ipk, then x,y are of 
    374377         ! size 1 (point/mooring data): this must be read onto the central grid point 
    375378         idvar  = iom_varid( sdjf%num, sdjf%clvar ) 
    376379         idmspc = iom_file ( sdjf%num )%ndims( idvar ) 
    377          IF( iom_file( sdjf%num )%luld( idvar ) )   idmspc = idmspc - 1 
    378          lmoor  = (  idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk  ) 
    379          ! 
    380          SELECT CASE( ipk ) 
    381          CASE(1) 
    382             IF( lk_c1d .AND. lmoor ) THEN 
    383                IF( sdjf%ln_tint ) THEN 
    384                   CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fdta(2,2,1,2), sdjf%nrec_a(1) ) 
    385                   CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,1,2),'Z',1. ) 
    386                ELSE 
    387                   CALL iom_get( sdjf%num, sdjf%clvar, sdjf%fnow(2,2,1  ), sdjf%nrec_a(1) ) 
    388                   CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,1  ),'Z',1. ) 
    389                ENDIF 
    390             ELSE 
    391                IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,1,2), sdjf%nrec_a(1) ) 
    392                ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,1  ), sdjf%nrec_a(1) ) 
    393                ENDIF 
    394             ENDIF 
    395          CASE DEFAULT 
    396             IF(lk_c1d .AND. lmoor ) THEN 
    397                IF( sdjf%ln_tint ) THEN 
    398                   CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) ) 
    399                   CALL lbc_lnk( 'fldread', sdjf%fdta(:,:,:,2),'Z',1. ) 
    400                ELSE 
    401                   CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fnow(2,2,:  ), sdjf%nrec_a(1) ) 
    402                   CALL lbc_lnk( 'fldread', sdjf%fnow(:,:,:  ),'Z',1. ) 
    403                ENDIF 
    404             ELSE 
    405                IF( sdjf%ln_tint ) THEN   ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fdta(:,:,:,2), sdjf%nrec_a(1) ) 
    406                ELSE                      ;   CALL iom_get( sdjf%num, ipdom, sdjf%clvar, sdjf%fnow(:,:,:  ), sdjf%nrec_a(1) ) 
    407                ENDIF 
    408             ENDIF 
    409          END SELECT 
    410       ENDIF 
    411       ! 
    412       sdjf%rotn(2) = .false.   ! vector not yet rotated 
     380         IF( iom_file( sdjf%num )%luld( idvar ) )   idmspc = idmspc - 1   ! id of the last spatial dimension 
     381         lmoor  = (  idmspc == 0 .OR. PRODUCT( iom_file( sdjf%num )%dimsz( 1:MAX(idmspc,1) ,idvar ) ) == ipk  )     
     382         ! 
     383         IF( lk_c1d .AND. lmoor ) THEN 
     384            CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, dta_alias(2,2,:), sdjf%nrec(1,iaa) )   ! jpdom_unknown -> no lbc_lnk 
     385            CALL lbc_lnk( 'fldread', dta_alias(:,:,:), 'T', 1., kfillmode = jpfillcopy ) 
     386         ELSE 
     387            CALL iom_get( sdjf%num,  jpdom_global, sdjf%clvar, dta_alias(:,:,:), sdjf%nrec(1,iaa),   & 
     388               &          sdjf%cltype, sdjf%zsgn, kfill = jpfillcopy ) 
     389         ENDIF 
     390      ENDIF 
     391      ! 
     392      sdjf%rotn(iaa) = .false.   ! vector not yet rotated 
    413393      ! 
    414394   END SUBROUTINE fld_get 
    … …  
    446426      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read_z  ! work space local data requiring vertical interpolation 
    447427      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   zdta_read_dz ! work space local data requiring vertical interpolation 
    448       CHARACTER(LEN=1),DIMENSION(3)             ::   clgrid 
     428      CHARACTER(LEN=1),DIMENSION(3)             ::   cltype 
    449429      LOGICAL                                   ::   lluld        ! is the variable using the unlimited dimension 
    450430      LOGICAL                                   ::   llzint       ! local value of ldzint 
    451431      !!--------------------------------------------------------------------- 
    452432      ! 
    453       clgrid = (/'t','u','v'/) 
     433      cltype = (/'t','u','v'/) 
    454434      ! 
    455435      ipi = SIZE( pdta, 1 ) 
    … …  
    486466         IF( ipkb /= ipk .OR. llzint ) THEN   ! boundary data not on model vertical grid : vertical interpolation 
    487467            ! 
    488             IF( ipk == jpk .AND. iom_varid(knum,'gdep'//clgrid(kgrd)) /= -1 .AND. iom_varid(knum,'e3'//clgrid(kgrd)) /= -1 ) THEN 
     468            IF( ipk == jpk .AND. iom_varid(knum,'gdep'//cltype(kgrd)) /= -1 .AND. iom_varid(knum,'e3'//cltype(kgrd)) /= -1 ) THEN 
    489469                
    490470               ALLOCATE( zdta_read(ipi,ipj,ipkb), zdta_read_z(ipi,ipj,ipkb), zdta_read_dz(ipi,ipj,ipkb) ) 
    491471                 
    492472               CALL fld_map_core( zz_read, kmap, zdta_read ) 
    493                CALL iom_get ( knum, jpdom_unknown, 'gdep'//clgrid(kgrd), zz_read )   ! read only once? Potential temporal evolution? 
     473               CALL iom_get ( knum, jpdom_unknown, 'gdep'//cltype(kgrd), zz_read )   ! read only once? Potential temporal evolution? 
    494474               CALL fld_map_core( zz_read, kmap, zdta_read_z ) 
    495                CALL iom_get ( knum, jpdom_unknown,   'e3'//clgrid(kgrd), zz_read )   ! read only once? Potential temporal evolution? 
     475               CALL iom_get ( knum, jpdom_unknown,   'e3'//cltype(kgrd), zz_read )   ! read only once? Potential temporal evolution? 
    496476               CALL fld_map_core( zz_read, kmap, zdta_read_dz ) 
    497477                
    … …  
    503483               IF( ipk /= jpk ) CALL ctl_stop( 'fld_map : this should be an impossible case...' ) 
    504484               WRITE(ctmp1,*) 'fld_map : vertical interpolation for bdy variable '//TRIM(cdvar)//' requires '  
    505                IF( iom_varid(knum, 'gdep'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//'gdep'//clgrid(kgrd)//' variable' ) 
    506                IF( iom_varid(knum,   'e3'//clgrid(kgrd)) == -1 ) CALL ctl_stop( ctmp1//  'e3'//clgrid(kgrd)//' variable' ) 
     485               IF( iom_varid(knum, 'gdep'//cltype(kgrd)) == -1 ) CALL ctl_stop( ctmp1//'gdep'//cltype(kgrd)//' variable' ) 
     486               IF( iom_varid(knum,   'e3'//cltype(kgrd)) == -1 ) CALL ctl_stop( ctmp1//  'e3'//cltype(kgrd)//' variable' ) 
    507487 
    508488            ENDIF 
    … …  
    632612               zdhalf(jk) = zdhalf(jk-1) + e3v(ji,jj,jk-1,Kmm) 
    633613               zdepth(jk) =          zcoef  * ( zdhalf(jk  ) + 0.5_wp * e3vw(ji,jj,jk,Kmm))  & 
    634                   &         + (1._wp-zcoef) * ( zdepth(jk-1) +          e3vw(ji,jj,jk,Kmm)) 
     614                     + (1._wp-zcoef) * ( zdepth(jk-1) +          e3vw(ji,jj,jk,Kmm)) 
    635615            END DO 
    636616         END SELECT 
    … …  
    727707      CHARACTER (LEN=100)          ::   clcomp       ! dummy weight name 
    728708      REAL(wp), DIMENSION(jpi,jpj) ::   utmp, vtmp   ! temporary arrays for vector rotation 
     709      REAL(wp), DIMENSION(:,:,:), POINTER ::   dta_u, dta_v    ! short cut 
    729710      !!--------------------------------------------------------------------- 
    730711      ! 
    … …  
    746727                  END DO 
    747728                  IF( iv > 0 ) THEN   ! fields ju and iv are two components which need to be rotated together 
     729                     IF( sd(ju)%ln_tint ) THEN   ;   dta_u => sd(ju)%fdta(:,:,:,jn)   ;   dta_v => sd(iv)%fdta(:,:,:,jn)  
     730                     ELSE                        ;   dta_u => sd(ju)%fnow(:,:,:   )   ;   dta_v => sd(iv)%fnow(:,:,:   ) 
     731                     ENDIF 
    748732                     DO jk = 1, SIZE( sd(ju)%fnow, 3 ) 
    749                         IF( sd(ju)%ln_tint )THEN 
    750                            CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->i', utmp(:,:) ) 
    751                            CALL rot_rep( sd(ju)%fdta(:,:,jk,jn), sd(iv)%fdta(:,:,jk,jn), 'T', 'en->j', vtmp(:,:) ) 
    752                            sd(ju)%fdta(:,:,jk,jn) = utmp(:,:)   ;   sd(iv)%fdta(:,:,jk,jn) = vtmp(:,:) 
    753                         ELSE  
    754                            CALL rot_rep( sd(ju)%fnow(:,:,jk   ), sd(iv)%fnow(:,:,jk   ), 'T', 'en->i', utmp(:,:) ) 
    755                            CALL rot_rep( sd(ju)%fnow(:,:,jk   ), sd(iv)%fnow(:,:,jk   ), 'T', 'en->j', vtmp(:,:) ) 
    756                            sd(ju)%fnow(:,:,jk   ) = utmp(:,:)   ;   sd(iv)%fnow(:,:,jk   ) = vtmp(:,:) 
    757                         ENDIF 
     733                        CALL rot_rep( dta_u(:,:,jk), dta_v(:,:,jk), 'T', 'en->i', utmp(:,:) ) 
     734                        CALL rot_rep( dta_u(:,:,jk), dta_v(:,:,jk), 'T', 'en->j', vtmp(:,:) ) 
     735                        dta_u(:,:,jk) = utmp(:,:)   ;   dta_v(:,:,jk) = vtmp(:,:) 
    758736                     END DO 
    759737                     sd(ju)%rotn(jn) = .TRUE.               ! vector was rotated  
    … …  
    801779 
    802780      ! current file parameters 
    803       IF( sdjf%cltype(1:4) == 'week' ) THEN          ! find the day of the beginning of the current week 
    804          isecwk = ksec_week( sdjf%cltype(6:8) )     ! seconds between the beginning of the week and half of current time step 
    805          llprevmt = isecwk > nsec_month               ! longer time since beginning of the current week than the current month 
     781      IF( sdjf%clftyp(1:4) == 'week' ) THEN         ! find the day of the beginning of the current week 
     782         isecwk = ksec_week( sdjf%clftyp(6:8) )     ! seconds between the beginning of the week and half of current time step 
     783         llprevmt = isecwk > nsec_month             ! longer time since beginning of the current week than the current month 
    806784         llprevyr = llprevmt .AND. nmonth == 1 
    807785         iyr = nyear  - COUNT((/llprevyr/)) 
    808786         imt = nmonth - COUNT((/llprevmt/)) + 12 * COUNT((/llprevyr/)) 
    809787         idy = nday + nmonth_len(nmonth-1) * COUNT((/llprevmt/)) - isecwk / idaysec 
    810          isecwk = nsec_year - isecwk              ! seconds between 00h jan 1st of current year and current week beginning 
     788         isecwk = nsec_year - isecwk                ! seconds between 00h jan 1st of current year and current week beginning 
    811789      ELSE 
    812790         iyr = nyear 
    … …  
    818796      ! previous file parameters 
    819797      IF( llprev ) THEN 
    820          IF( sdjf%cltype(1:4) == 'week'    ) THEN     ! find the day of the beginning of previous week 
    821             isecwk = isecwk + 7 * idaysec         ! seconds between the beginning of previous week and half of the time step 
    822             llprevmt = isecwk > nsec_month            ! longer time since beginning of the previous week than the current month 
     798         IF( sdjf%clftyp(1:4) == 'week'    ) THEN   ! find the day of the beginning of previous week 
     799            isecwk = isecwk + 7 * idaysec           ! seconds between the beginning of previous week and half of the time step 
     800            llprevmt = isecwk > nsec_month          ! longer time since beginning of the previous week than the current month 
    823801            llprevyr = llprevmt .AND. nmonth == 1 
    824802            iyr = nyear  - COUNT((/llprevyr/)) 
    825803            imt = nmonth - COUNT((/llprevmt/)) + 12 * COUNT((/llprevyr/)) 
    826804            idy = nday + nmonth_len(nmonth-1) * COUNT((/llprevmt/)) - isecwk / idaysec 
    827             isecwk = nsec_year - isecwk           ! seconds between 00h jan 1st of current year and previous week beginning 
     805            isecwk = nsec_year - isecwk             ! seconds between 00h jan 1st of current year and previous week beginning 
    828806         ELSE 
    829             idy = nday   - COUNT((/ sdjf%cltype == 'daily'                 /)) 
    830             imt = nmonth - COUNT((/ sdjf%cltype == 'monthly' .OR. idy == 0 /)) 
    831             iyr = nyear  - COUNT((/ sdjf%cltype == 'yearly'  .OR. imt == 0 /)) 
     807            idy = nday   - COUNT((/ sdjf%clftyp == 'daily'                 /)) 
     808            imt = nmonth - COUNT((/ sdjf%clftyp == 'monthly' .OR. idy == 0 /)) 
     809            iyr = nyear  - COUNT((/ sdjf%clftyp == 'yearly'  .OR. imt == 0 /)) 
    832810            IF( idy == 0 ) idy = nmonth_len(imt) 
    833811            IF( imt == 0 ) imt = 12 
    … …  
    838816      ! next file parameters 
    839817      IF( llnext ) THEN 
    840          IF( sdjf%cltype(1:4) == 'week'    ) THEN     ! find the day of the beginning of next week 
    841             isecwk = 7 * idaysec - isecwk         ! seconds between half of the time step and the beginning of next week 
     818         IF( sdjf%clftyp(1:4) == 'week'    ) THEN   ! find the day of the beginning of next week 
     819            isecwk = 7 * idaysec - isecwk           ! seconds between half of the time step and the beginning of next week 
    842820            llnextmt = isecwk > ( nmonth_len(nmonth)*idaysec - nsec_month )   ! larger than the seconds to the end of the month 
    843821            llnextyr = llnextmt .AND. nmonth == 12 
    … …  
    845823            imt = nmonth + COUNT((/llnextmt/)) - 12 * COUNT((/llnextyr/)) 
    846824            idy = nday - nmonth_len(nmonth) * COUNT((/llnextmt/)) + isecwk / idaysec + 1 
    847             isecwk = nsec_year + isecwk           ! seconds between 00h jan 1st of current year and next week beginning 
     825            isecwk = nsec_year + isecwk             ! seconds between 00h jan 1st of current year and next week beginning 
    848826         ELSE 
    849             idy = nday   + COUNT((/ sdjf%cltype == 'daily'                                 /)) 
    850             imt = nmonth + COUNT((/ sdjf%cltype == 'monthly' .OR. idy > nmonth_len(nmonth) /)) 
    851             iyr = nyear  + COUNT((/ sdjf%cltype == 'yearly'  .OR. imt == 13                /)) 
     827            idy = nday   + COUNT((/ sdjf%clftyp == 'daily'                                 /)) 
     828            imt = nmonth + COUNT((/ sdjf%clftyp == 'monthly' .OR. idy > nmonth_len(nmonth) /)) 
     829            iyr = nyear  + COUNT((/ sdjf%clftyp == 'yearly'  .OR. imt == 13                /)) 
    852830            IF( idy > nmonth_len(nmonth) )   idy = 1 
    853831            IF( imt == 13                )   imt = 1 
    … …  
    866844      IF    ( NINT(sdjf%freqh) == -12 ) THEN            ;   ireclast = 1    ! yearly mean: consider only 1 record 
    867845      ELSEIF( NINT(sdjf%freqh) ==  -1 ) THEN                                ! monthly mean: 
    868          IF(     sdjf%cltype      == 'monthly' ) THEN   ;   ireclast = 1    !  consider that the file has  1 record 
     846         IF(     sdjf%clftyp      == 'monthly' ) THEN   ;   ireclast = 1    !  consider that the file has  1 record 
    869847         ELSE                                           ;   ireclast = 12   !  consider that the file has 12 record 
    870848         ENDIF 
    871849      ELSE                                                                  ! higher frequency mean (in hours) 
    872          IF(     sdjf%cltype      == 'monthly' ) THEN   ;   ireclast = NINT( 24. * REAL(nmonth_len(indexmt), wp) / sdjf%freqh ) 
    873          ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   ireclast = NINT( 24. * 7.                            / sdjf%freqh ) 
    874          ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   ireclast = NINT( 24.                                 / sdjf%freqh ) 
     850         IF(     sdjf%clftyp      == 'monthly' ) THEN   ;   ireclast = NINT( 24. * REAL(nmonth_len(indexmt), wp) / sdjf%freqh ) 
     851         ELSEIF( sdjf%clftyp(1:4) == 'week'    ) THEN   ;   ireclast = NINT( 24. * 7.                            / sdjf%freqh ) 
     852         ELSEIF( sdjf%clftyp      == 'daily'   ) THEN   ;   ireclast = NINT( 24.                                 / sdjf%freqh ) 
    875853         ELSE                                           ;   ireclast = NINT( 24. * REAL( nyear_len(indexyr), wp) / sdjf%freqh ) 
    876854         ENDIF 
    … …  
    890868         sdjf%nrecsec(1) = sdjf%nrecsec(0) + nyear_len( indexyr ) * idaysec 
    891869      ELSEIF( NINT(sdjf%freqh) ==  -1 ) THEN                                     ! monthly mean: 
    892          IF(     sdjf%cltype      == 'monthly' ) THEN                            !    monthly file 
     870         IF(     sdjf%clftyp      == 'monthly' ) THEN                            !    monthly file 
    893871            sdjf%nrecsec(0   ) = nsec1jan000 + nmonth_beg(indexmt  ) 
    894872            sdjf%nrecsec(1   ) = nsec1jan000 + nmonth_beg(indexmt+1) 
    … …  
    898876         ENDIF 
    899877      ELSE                                                                       ! higher frequency mean (in hours) 
    900          IF(     sdjf%cltype      == 'monthly' ) THEN   ;   istart = nsec1jan000 + nmonth_beg(indexmt) 
    901          ELSEIF( sdjf%cltype(1:4) == 'week'    ) THEN   ;   istart = nsec1jan000 + isecwk 
    902          ELSEIF( sdjf%cltype      == 'daily'   ) THEN   ;   istart = nsec1jan000 + nmonth_beg(indexmt) + ( idy - 1 ) * idaysec 
     878         IF(     sdjf%clftyp      == 'monthly' ) THEN   ;   istart = nsec1jan000 + nmonth_beg(indexmt) 
     879         ELSEIF( sdjf%clftyp(1:4) == 'week'    ) THEN   ;   istart = nsec1jan000 + isecwk 
     880         ELSEIF( sdjf%clftyp      == 'daily'   ) THEN   ;   istart = nsec1jan000 + nmonth_beg(indexmt) + ( idy - 1 ) * idaysec 
    903881         ELSEIF( indexyr          == 0         ) THEN   ;   istart = nsec1jan000 - nyear_len( 0 ) * idaysec 
    904882         ELSEIF( indexyr          == 2         ) THEN   ;   istart = nsec1jan000 + nyear_len( 1 ) * idaysec 
    … …  
    941919      IF( sdjf%num <= 0 .OR. .NOT. sdjf%ln_clim  ) THEN 
    942920         IF( sdjf%num > 0 )   CALL iom_close( sdjf%num )   ! close file if already open 
    943          CALL iom_open( sdjf%clname, sdjf%num, ldstop = llstop, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 ) 
     921         CALL iom_open( sdjf%clname, sdjf%num, ldstop = llstop, ldiof = LEN_TRIM(sdjf%wgtname) > 0 ) 
    944922      ENDIF 
    945923      ! 
    … …  
    963941         ENDIF 
    964942         ! 
    965          CALL iom_open( sdjf%clname, sdjf%num, ldiof = LEN(TRIM(sdjf%wgtname)) > 0 )    
     943         CALL iom_open( sdjf%clname, sdjf%num, ldiof = LEN_TRIM(sdjf%wgtname) > 0 )    
    966944         ! 
    967945      ENDIF 
    … …  
    996974         sdf(jf)%ln_tint    = sdf_n(jf)%ln_tint 
    997975         sdf(jf)%ln_clim    = sdf_n(jf)%ln_clim 
    998          sdf(jf)%cltype     = sdf_n(jf)%cltype 
     976         sdf(jf)%clftyp     = sdf_n(jf)%clftyp 
     977         sdf(jf)%cltype     = 'T'   ! by default don't do any call to lbc_lnk in iom_get 
     978         sdf(jf)%zsgn       = 1.    ! by default don't do change signe across the north fold 
    999979         sdf(jf)%num        = -1 
     980         sdf(jf)%nbb        = 1  ! start with before data in 1 
     981         sdf(jf)%naa        = 2  ! start with after  data in 2 
    1000982         sdf(jf)%wgtname    = " " 
    1001983         IF( LEN( TRIM(sdf_n(jf)%wname) ) > 0 )   sdf(jf)%wgtname = TRIM( cdir )//sdf_n(jf)%wname 
    … …  
    1004986         sdf(jf)%vcomp      = sdf_n(jf)%vcomp 
    1005987         sdf(jf)%rotn(:)    = .TRUE.   ! pretend to be rotated -> won't try to rotate data before the first call to fld_get 
    1006          IF( sdf(jf)%cltype(1:4) == 'week' .AND. nn_leapy == 0  )   & 
     988         IF( sdf(jf)%clftyp(1:4) == 'week' .AND. nn_leapy == 0  )   & 
    1007989            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs nn_leapy = 1') 
    1008          IF( sdf(jf)%cltype(1:4) == 'week' .AND. sdf(jf)%ln_clim )   & 
     990         IF( sdf(jf)%clftyp(1:4) == 'week' .AND. sdf(jf)%ln_clim )   & 
    1009991            &   CALL ctl_stop('fld_clopn: weekly file ('//TRIM(sdf(jf)%clrootname)//') needs ln_clim = .FALSE.') 
    1010992         sdf(jf)%nreclast   = -1 ! Set to non zero default value to avoid errors, is updated to meaningful value during fld_clopn 
    … …  
    10321014            WRITE(numout,*) '         weights: '        , TRIM( sdf(jf)%wgtname    ),   & 
    10331015               &                  '   pairing: '        , TRIM( sdf(jf)%vcomp      ),   & 
    1034                &                  '   data type: '      ,       sdf(jf)%cltype      ,   & 
     1016               &                  '   data type: '      ,       sdf(jf)%clftyp      ,   & 
    10351017               &                  '   land/sea mask:'   , TRIM( sdf(jf)%lsmname    ) 
    10361018            call flush(numout) 
    … …  
    10501032      !!---------------------------------------------------------------------- 
    10511033      TYPE( FLD ), INTENT(in   ) ::   sd        ! field with name of weights file 
    1052       INTEGER    , INTENT(inout) ::   kwgt      ! index of weights 
     1034      INTEGER    , INTENT(  out) ::   kwgt      ! index of weights 
    10531035      ! 
    10541036      INTEGER ::   kw, nestid   ! local integer 
    1055       LOGICAL ::   found        ! local logical 
    10561037      !!---------------------------------------------------------------------- 
    10571038      ! 
    10581039      !! search down linked list  
    10591040      !! weights filename is either present or we hit the end of the list 
    1060       found = .FALSE. 
    10611041      ! 
    10621042      !! because agrif nest part of filenames are now added in iom_open 
    … …  
    10681048#endif 
    10691049      DO kw = 1, nxt_wgt-1 
    1070          IF( TRIM(ref_wgts(kw)%wgtname) == TRIM(sd%wgtname) .AND. & 
    1071              ref_wgts(kw)%nestid == nestid) THEN 
     1050         IF( ref_wgts(kw)%wgtname == sd%wgtname .AND. & 
     1051             ref_wgts(kw)%nestid  == nestid) THEN 
    10721052            kwgt = kw 
    1073             found = .TRUE. 
    1074             EXIT 
     1053            RETURN 
    10751054         ENDIF 
    10761055      END DO 
    1077       IF( .NOT.found ) THEN 
    1078          kwgt = nxt_wgt 
    1079          CALL fld_weight( sd ) 
    1080       ENDIF 
     1056      kwgt = nxt_wgt 
     1057      CALL fld_weight( sd ) 
    10811058      ! 
    10821059   END SUBROUTINE wgt_list 
    … …  
    11211098      TYPE( FLD ), INTENT(in) ::   sd   ! field with name of weights file 
    11221099      !! 
    1123       INTEGER ::   jn         ! dummy loop indices 
     1100      INTEGER ::   ji,jj,jn   ! dummy loop indices 
    11241101      INTEGER ::   inum       ! local logical unit 
    11251102      INTEGER ::   id         ! local variable id 
    … …  
    11271104      INTEGER ::   zwrap      ! local integer 
    11281105      LOGICAL ::   cyclical   !  
    1129       CHARACTER (len=5) ::   aname   ! 
    1130       INTEGER , DIMENSION(:), ALLOCATABLE ::   ddims 
    1131       INTEGER,  DIMENSION(jpi,jpj) ::   data_src 
     1106      CHARACTER (len=5) ::   clname   ! 
     1107      INTEGER , DIMENSION(4) ::   ddims 
     1108      INTEGER                ::   isrc 
    11321109      REAL(wp), DIMENSION(jpi,jpj) ::   data_tmp 
    11331110      !!---------------------------------------------------------------------- 
    … …  
    11421119      !! current weights file 
    11431120 
    1144       !! open input data file (non-model grid) 
    1145       CALL iom_open( sd%clname, inum, ldiof =  LEN(TRIM(sd%wgtname)) > 0 ) 
    1146  
    1147       !! get dimensions: we consider 2D data as 3D data with vertical dim size = 1 
    1148       IF( SIZE(sd%fnow, 3) > 0 ) THEN 
    1149          ALLOCATE( ddims(4) ) 
    1150       ELSE 
    1151          ALLOCATE( ddims(3) ) 
    1152       ENDIF 
    1153       id = iom_varid( inum, sd%clvar, ddims ) 
    1154  
    1155       !! close it 
    1156       CALL iom_close( inum ) 
     1121      !! get data grid dimensions 
     1122      id = iom_varid( sd%num, sd%clvar, ddims ) 
    11571123 
    11581124      !! now open the weights file 
    1159  
    11601125      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights 
    11611126      IF( inum > 0 ) THEN 
    … …  
    11931158         !! two possible cases: bilinear (4 weights) or bicubic (16 weights) 
    11941159         id = iom_varid(inum, 'src05', ldstop=.FALSE.) 
    1195          IF( id <= 0) THEN 
    1196             ref_wgts(nxt_wgt)%numwgt = 4 
    1197          ELSE 
    1198             ref_wgts(nxt_wgt)%numwgt = 16 
    1199          ENDIF 
    1200  
    1201          ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(jpi,jpj,4) ) 
    1202          ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(jpi,jpj,4) ) 
    1203          ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(jpi,jpj,ref_wgts(nxt_wgt)%numwgt) ) 
     1160         IF( id <= 0 ) THEN   ;   ref_wgts(nxt_wgt)%numwgt = 4 
     1161         ELSE                 ;   ref_wgts(nxt_wgt)%numwgt = 16 
     1162         ENDIF 
     1163 
     1164         ALLOCATE( ref_wgts(nxt_wgt)%data_jpi(Nis0:Nie0,Njs0:Nje0,4) ) 
     1165         ALLOCATE( ref_wgts(nxt_wgt)%data_jpj(Nis0:Nie0,Njs0:Nje0,4) ) 
     1166         ALLOCATE( ref_wgts(nxt_wgt)%data_wgt(Nis0:Nie0,Njs0:Nje0,ref_wgts(nxt_wgt)%numwgt) ) 
    12041167 
    12051168         DO jn = 1,4 
    1206             aname = ' ' 
    1207             WRITE(aname,'(a3,i2.2)') 'src',jn 
    1208             data_tmp(:,:) = 0 
    1209             CALL iom_get ( inum, jpdom_data, aname, data_tmp(:,:) ) 
    1210             data_src(:,:) = INT(data_tmp(:,:)) 
    1211             ref_wgts(nxt_wgt)%data_jpj(:,:,jn) = 1 + (data_src(:,:)-1) / ref_wgts(nxt_wgt)%ddims(1) 
    1212             ref_wgts(nxt_wgt)%data_jpi(:,:,jn) = data_src(:,:) - ref_wgts(nxt_wgt)%ddims(1)*(ref_wgts(nxt_wgt)%data_jpj(:,:,jn)-1) 
     1169            WRITE(clname,'(a3,i2.2)') 'src',jn 
     1170            CALL iom_get ( inum, jpdom_global, clname, data_tmp(:,:), cd_type = 'Z' )   !  no call to lbc_lnk 
     1171            DO_2D( 0, 0, 0, 0 ) 
     1172               isrc = NINT(data_tmp(ji,jj)) - 1 
     1173               ref_wgts(nxt_wgt)%data_jpi(ji,jj,jn) = 1 + MOD(isrc,  ref_wgts(nxt_wgt)%ddims(1)) 
     1174               ref_wgts(nxt_wgt)%data_jpj(ji,jj,jn) = 1 +     isrc / ref_wgts(nxt_wgt)%ddims(1) 
     1175            END_2D 
    12131176         END DO 
    12141177 
    12151178         DO jn = 1, ref_wgts(nxt_wgt)%numwgt 
    1216             aname = ' ' 
    1217             WRITE(aname,'(a3,i2.2)') 'wgt',jn 
    1218             ref_wgts(nxt_wgt)%data_wgt(:,:,jn) = 0.0 
    1219             CALL iom_get ( inum, jpdom_data, aname, ref_wgts(nxt_wgt)%data_wgt(:,:,jn) ) 
     1179            WRITE(clname,'(a3,i2.2)') 'wgt',jn 
     1180            CALL iom_get ( inum, jpdom_global, clname, data_tmp(:,:), cd_type = 'Z' )   !  no call to lbc_lnk 
     1181            DO_2D( 0, 0, 0, 0 ) 
     1182               ref_wgts(nxt_wgt)%data_wgt(ji,jj,jn) = data_tmp(ji,jj) 
     1183            END_2D 
    12201184         END DO 
    12211185         CALL iom_close (inum) 
    12221186  
    12231187         ! find min and max indices in grid 
    1224          ref_wgts(nxt_wgt)%botleft(1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:)) 
    1225          ref_wgts(nxt_wgt)%botleft(2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:)) 
     1188         ref_wgts(nxt_wgt)%botleft( 1) = MINVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:)) 
     1189         ref_wgts(nxt_wgt)%botleft( 2) = MINVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:)) 
    12261190         ref_wgts(nxt_wgt)%topright(1) = MAXVAL(ref_wgts(nxt_wgt)%data_jpi(:,:,:)) 
    12271191         ref_wgts(nxt_wgt)%topright(2) = MAXVAL(ref_wgts(nxt_wgt)%data_jpj(:,:,:)) 
    … …  
    12471211         CALL ctl_stop( '    fld_weight : unable to read the file ' ) 
    12481212      ENDIF 
    1249  
    1250       DEALLOCATE (ddims ) 
    12511213      ! 
    12521214   END SUBROUTINE fld_weight 
    … …  
    12811243      SELECT CASE( SIZE(zfieldo(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),3) ) 
    12821244      CASE(1) 
    1283          CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1), 1, rec1_lsm, recn_lsm) 
     1245         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1),   & 
     1246            &          1, kstart = rec1_lsm, kcount = recn_lsm) 
    12841247      CASE DEFAULT 
    1285          CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:), 1, rec1_lsm, recn_lsm) 
     1248         CALL iom_get( inum, jpdom_unknown, 'LSM', zslmec1(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),   & 
     1249            &          1, kstart = rec1_lsm, kcount = recn_lsm) 
    12861250      END SELECT 
    12871251      CALL iom_close( inum ) 
    … …  
    13261290      !!      D. Delrosso INGV 
    13271291      !!----------------------------------------------------------------------  
    1328       INTEGER                      , INTENT(in   ) :: ileni,ilenj   ! lengths  
    1329       REAL, DIMENSION (ileni,ilenj), INTENT(in   ) :: zfieldn       ! array of forcing field with undeff for land points 
    1330       REAL, DIMENSION (ileni,ilenj), INTENT(  out) :: zfield        ! array of forcing field 
    1331       ! 
    1332       REAL  , DIMENSION (ileni,ilenj)   :: zmat1, zmat2, zmat3, zmat4  ! local arrays  
    1333       REAL  , DIMENSION (ileni,ilenj)   :: zmat5, zmat6, zmat7, zmat8  !   -     -  
    1334       REAL  , DIMENSION (ileni,ilenj)   :: zlsm2d                      !   -     -  
    1335       REAL  , DIMENSION (ileni,ilenj,8) :: zlsm3d                      !   -     - 
    1336       LOGICAL, DIMENSION (ileni,ilenj,8) :: ll_msknan3d                 ! logical mask for undeff detection 
    1337       LOGICAL, DIMENSION (ileni,ilenj)   :: ll_msknan2d                 ! logical mask for undeff detection 
     1292      INTEGER                          , INTENT(in   ) :: ileni,ilenj   ! lengths  
     1293      REAL(wp), DIMENSION (ileni,ilenj), INTENT(in   ) :: zfieldn       ! array of forcing field with undeff for land points 
     1294      REAL(wp), DIMENSION (ileni,ilenj), INTENT(  out) :: zfield        ! array of forcing field 
     1295      ! 
     1296      REAL(wp) , DIMENSION (ileni,ilenj)   :: zmat1, zmat2, zmat3, zmat4  ! local arrays  
     1297      REAL(wp) , DIMENSION (ileni,ilenj)   :: zmat5, zmat6, zmat7, zmat8  !   -     -  
     1298      REAL(wp) , DIMENSION (ileni,ilenj)   :: zlsm2d                      !   -     -  
     1299      REAL(wp) , DIMENSION (ileni,ilenj,8) :: zlsm3d                      !   -     - 
     1300      LOGICAL  , DIMENSION (ileni,ilenj,8) :: ll_msknan3d                 ! logical mask for undeff detection 
     1301      LOGICAL  , DIMENSION (ileni,ilenj)   :: ll_msknan2d                 ! logical mask for undeff detection 
    13381302      !!----------------------------------------------------------------------  
    13391303      zmat8 = eoshift( zfieldn , SHIFT=-1 , BOUNDARY = (/zfieldn(:,1)/)     , DIM=2 ) 
    … …  
    13561320 
    13571321 
    1358    SUBROUTINE fld_interp( num, clvar, kw, kk, dta,  & 
    1359                           &         nrec, lsmfile)       
     1322   SUBROUTINE fld_interp( num, clvar, kw, kk, dta, nrec, lsmfile)       
    13601323      !!--------------------------------------------------------------------- 
    13611324      !!                    ***  ROUTINE fld_interp  *** 
    … …  
    13751338      INTEGER, DIMENSION(3) ::   rec1_lsm, recn_lsm   ! temporary arrays for start and length in case of seaoverland 
    13761339      INTEGER ::   ii_lsm1,ii_lsm2,ij_lsm1,ij_lsm2    ! temporary indices 
    1377       INTEGER ::   jk, jn, jm, jir, jjr               ! loop counters 
     1340      INTEGER ::   ji, jj, jk, jn, jir, jjr           ! loop counters 
     1341      INTEGER ::   ipk 
    13781342      INTEGER ::   ni, nj                             ! lengths 
    13791343      INTEGER ::   jpimin,jpiwid                      ! temporary indices 
    … …  
    13861350      REAL(wp),DIMENSION(:,:,:), ALLOCATABLE ::   ztmp_fly_dta                 ! local array of values on input grid      
    13871351      !!---------------------------------------------------------------------- 
     1352      ipk = SIZE(dta, 3) 
    13881353      ! 
    13891354      !! for weighted interpolation we have weights at four corners of a box surrounding  
    … …  
    14151380 
    14161381 
    1417       IF( LEN( TRIM(lsmfile) ) > 0 ) THEN 
     1382      IF( LEN_TRIM(lsmfile) > 0 ) THEN 
    14181383      !! indeces for ztmp_fly_dta 
    14191384      ! -------------------------- 
    … …  
    14451410         CASE(1) 
    14461411              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,1),   & 
    1447                  &                                                                nrec, rec1_lsm, recn_lsm) 
     1412                 &          nrec, kstart = rec1_lsm, kcount = recn_lsm) 
    14481413         CASE DEFAULT 
    14491414              CALL iom_get( num, jpdom_unknown, clvar, ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),   & 
    1450                  &                                                                nrec, rec1_lsm, recn_lsm) 
     1415                 &          nrec, kstart = rec1_lsm, kcount = recn_lsm) 
    14511416         END SELECT 
    14521417         CALL apply_seaoverland(lsmfile,ztmp_fly_dta(jpi1_lsm:jpi2_lsm,jpj1_lsm:jpj2_lsm,:),                  & 
    … …  
    14681433          
    14691434         ref_wgts(kw)%fly_dta(:,:,:) = 0.0 
    1470          CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, rec1, recn) 
     1435         CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%fly_dta(jpi1:jpi2,jpj1:jpj2,:), nrec, kstart = rec1, kcount = recn) 
    14711436      ENDIF 
    14721437       
    … …  
    14741439      !! first four weights common to both bilinear and bicubic 
    14751440      !! data_jpi, data_jpj have already been shifted to (1,1) corresponding to botleft 
    1476       !! note that we have to offset by 1 into fly_dta array because of halo 
    1477       dta(:,:,:) = 0.0 
    1478       DO jk = 1,4 
    1479         DO jn = 1, jpj 
    1480           DO jm = 1,jpi 
    1481             ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    1482             nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
    1483             dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk) * ref_wgts(kw)%fly_dta(ni+1,nj+1,:) 
    1484           END DO 
    1485         END DO 
     1441      !! note that we have to offset by 1 into fly_dta array because of halo added to fly_dta (rec1 definition) 
     1442      dta(:,:,:) = 0._wp 
     1443      DO jn = 1,4 
     1444         DO_3D( 0, 0, 0, 0, 1,ipk ) 
     1445            ni = ref_wgts(kw)%data_jpi(ji,jj,jn) + 1 
     1446            nj = ref_wgts(kw)%data_jpj(ji,jj,jn) + 1 
     1447            dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn) * ref_wgts(kw)%fly_dta(ni,nj,jk) 
     1448         END_3D 
    14861449      END DO 
    14871450 
    14881451      IF(ref_wgts(kw)%numwgt .EQ. 16) THEN 
    14891452 
    1490         !! fix up halo points that we couldnt read from file 
    1491         IF( jpi1 == 2 ) THEN 
    1492            ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:) 
    1493         ENDIF 
    1494         IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 
    1495            ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:) 
    1496         ENDIF 
    1497         IF( jpj1 == 2 ) THEN 
    1498            ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:) 
    1499         ENDIF 
    1500         IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .lt. jpjwid+2 ) THEN 
    1501            ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:) 
    1502         ENDIF 
    1503  
    1504         !! if data grid is cyclic we can do better on east-west edges 
    1505         !! but have to allow for whether first and last columns are coincident 
    1506         IF( ref_wgts(kw)%cyclic ) THEN 
    1507            rec1(2) = MAX( jpjmin-1, 1 ) 
    1508            recn(1) = 1 
    1509            recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 ) 
    1510            jpj1 = 2 + rec1(2) - jpjmin 
    1511            jpj2 = jpj1 + recn(2) - 1 
    1512            IF( jpi1 == 2 ) THEN 
    1513               rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap 
    1514               CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 
    1515               ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 
    1516            ENDIF 
    1517            IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 
    1518               rec1(1) = 1 + ref_wgts(kw)%overlap 
    1519               CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, rec1, recn) 
    1520               ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 
    1521            ENDIF 
    1522         ENDIF 
    1523  
    1524         ! gradient in the i direction 
    1525         DO jk = 1,4 
    1526           DO jn = 1, jpj 
    1527             DO jm = 1,jpi 
    1528               ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    1529               nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
    1530               dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+4) * 0.5 *         & 
    1531                                (ref_wgts(kw)%fly_dta(ni+2,nj+1,:) - ref_wgts(kw)%fly_dta(ni,nj+1,:)) 
    1532             END DO 
    1533           END DO 
    1534         END DO 
    1535  
    1536         ! gradient in the j direction 
    1537         DO jk = 1,4 
    1538           DO jn = 1, jpj 
    1539             DO jm = 1,jpi 
    1540               ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    1541               nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
    1542               dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+8) * 0.5 *         & 
    1543                                (ref_wgts(kw)%fly_dta(ni+1,nj+2,:) - ref_wgts(kw)%fly_dta(ni+1,nj,:)) 
    1544             END DO 
    1545           END DO 
    1546         END DO 
    1547  
    1548          ! gradient in the ij direction 
    1549          DO jk = 1,4 
    1550             DO jn = 1, jpj 
    1551                DO jm = 1,jpi 
    1552                   ni = ref_wgts(kw)%data_jpi(jm,jn,jk) 
    1553                   nj = ref_wgts(kw)%data_jpj(jm,jn,jk) 
    1554                   dta(jm,jn,:) = dta(jm,jn,:) + ref_wgts(kw)%data_wgt(jm,jn,jk+12) * 0.25 * ( & 
    1555                                (ref_wgts(kw)%fly_dta(ni+2,nj+2,:) - ref_wgts(kw)%fly_dta(ni  ,nj+2,:)) -   & 
    1556                                (ref_wgts(kw)%fly_dta(ni+2,nj  ,:) - ref_wgts(kw)%fly_dta(ni  ,nj  ,:))) 
    1557                END DO 
    1558             END DO 
     1453         !! fix up halo points that we couldnt read from file 
     1454         IF( jpi1 == 2 ) THEN 
     1455            ref_wgts(kw)%fly_dta(jpi1-1,:,:) = ref_wgts(kw)%fly_dta(jpi1,:,:) 
     1456         ENDIF 
     1457         IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 
     1458            ref_wgts(kw)%fly_dta(jpi2+1,:,:) = ref_wgts(kw)%fly_dta(jpi2,:,:) 
     1459         ENDIF 
     1460         IF( jpj1 == 2 ) THEN 
     1461            ref_wgts(kw)%fly_dta(:,jpj1-1,:) = ref_wgts(kw)%fly_dta(:,jpj1,:) 
     1462         ENDIF 
     1463         IF( jpj2 + jpjmin - 1 == ref_wgts(kw)%ddims(2)+1 .AND. jpj2 .LT. jpjwid+2 ) THEN 
     1464            ref_wgts(kw)%fly_dta(:,jpj2+1,:) = 2.0*ref_wgts(kw)%fly_dta(:,jpj2,:) - ref_wgts(kw)%fly_dta(:,jpj2-1,:) 
     1465         ENDIF 
     1466          
     1467         !! if data grid is cyclic we can do better on east-west edges 
     1468         !! but have to allow for whether first and last columns are coincident 
     1469         IF( ref_wgts(kw)%cyclic ) THEN 
     1470            rec1(2) = MAX( jpjmin-1, 1 ) 
     1471            recn(1) = 1 
     1472            recn(2) = MIN( jpjwid+2, ref_wgts(kw)%ddims(2)-rec1(2)+1 ) 
     1473            jpj1 = 2 + rec1(2) - jpjmin 
     1474            jpj2 = jpj1 + recn(2) - 1 
     1475            IF( jpi1 == 2 ) THEN 
     1476               rec1(1) = ref_wgts(kw)%ddims(1) - ref_wgts(kw)%overlap 
     1477               CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, kstart = rec1, kcount = recn) 
     1478               ref_wgts(kw)%fly_dta(jpi1-1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 
     1479            ENDIF 
     1480            IF( jpi2 + jpimin - 1 == ref_wgts(kw)%ddims(1)+1 ) THEN 
     1481               rec1(1) = 1 + ref_wgts(kw)%overlap 
     1482               CALL iom_get( num, jpdom_unknown, clvar, ref_wgts(kw)%col(:,jpj1:jpj2,:), nrec, kstart = rec1, kcount = recn) 
     1483               ref_wgts(kw)%fly_dta(jpi2+1,jpj1:jpj2,:) = ref_wgts(kw)%col(1,jpj1:jpj2,:) 
     1484            ENDIF 
     1485         ENDIF 
     1486         ! 
     1487!!$         DO jn = 1,4 
     1488!!$            DO_3D( 0, 0, 0, 0, 1,ipk ) 
     1489!!$               ni = ref_wgts(kw)%data_jpi(ji,jj,jn) + 1 
     1490!!$               nj = ref_wgts(kw)%data_jpj(ji,jj,jn) + 1 
     1491!!$               dta(ji,jj,jk) = dta(ji,jj,jk)   & 
     1492!!$                  ! gradient in the i direction 
     1493!!$                  &            + ref_wgts(kw)%data_wgt(ji,jj,jn+4) * 0.5_wp *                                    & 
     1494!!$                  &                (ref_wgts(kw)%fly_dta(ni+1,nj  ,jk) - ref_wgts(kw)%fly_dta(ni-1,nj  ,jk))     & 
     1495!!$                  ! gradient in the j direction 
     1496!!$                  &            + ref_wgts(kw)%data_wgt(ji,jj,jn+8) * 0.5_wp *                                    & 
     1497!!$                  &                (ref_wgts(kw)%fly_dta(ni  ,nj+1,jk) - ref_wgts(kw)%fly_dta(ni  ,nj-1,jk))     & 
     1498!!$                  ! gradient in the ij direction 
     1499!!$                  &            + ref_wgts(kw)%data_wgt(ji,jj,jn+12) * 0.25_wp *                                  & 
     1500!!$                  &               ((ref_wgts(kw)%fly_dta(ni+1,nj+1,jk) - ref_wgts(kw)%fly_dta(ni-1,nj+1,jk)) -   & 
     1501!!$                  &                (ref_wgts(kw)%fly_dta(ni+1,nj-1,jk) - ref_wgts(kw)%fly_dta(ni-1,nj-1,jk))) 
     1502!!$            END_3D 
     1503!!$         END DO 
     1504         ! 
     1505         DO jn = 1,4 
     1506            DO_3D( 0, 0, 0, 0, 1,ipk ) 
     1507               ni = ref_wgts(kw)%data_jpi(ji,jj,jn) 
     1508               nj = ref_wgts(kw)%data_jpj(ji,jj,jn) 
     1509               ! gradient in the i direction 
     1510               dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn+4) * 0.5_wp *         & 
     1511                  &                (ref_wgts(kw)%fly_dta(ni+2,nj+1,jk) - ref_wgts(kw)%fly_dta(ni  ,nj+1,jk)) 
     1512            END_3D 
     1513         END DO 
     1514         DO jn = 1,4 
     1515            DO_3D( 0, 0, 0, 0, 1,ipk ) 
     1516               ni = ref_wgts(kw)%data_jpi(ji,jj,jn) 
     1517               nj = ref_wgts(kw)%data_jpj(ji,jj,jn) 
     1518               ! gradient in the j direction 
     1519               dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn+8) * 0.5_wp *         & 
     1520                  &                (ref_wgts(kw)%fly_dta(ni+1,nj+2,jk) - ref_wgts(kw)%fly_dta(ni+1,nj  ,jk)) 
     1521            END_3D 
     1522         END DO 
     1523         DO jn = 1,4 
     1524            DO_3D( 0, 0, 0, 0, 1,ipk ) 
     1525               ni = ref_wgts(kw)%data_jpi(ji,jj,jn) 
     1526               nj = ref_wgts(kw)%data_jpj(ji,jj,jn) 
     1527               ! gradient in the ij direction 
     1528               dta(ji,jj,jk) = dta(ji,jj,jk) + ref_wgts(kw)%data_wgt(ji,jj,jn+12) * 0.25_wp * (     & 
     1529                  &                (ref_wgts(kw)%fly_dta(ni+2,nj+2,jk) - ref_wgts(kw)%fly_dta(ni  ,nj+2,jk)) -   & 
     1530                  &                (ref_wgts(kw)%fly_dta(ni+2,nj  ,jk) - ref_wgts(kw)%fly_dta(ni  ,nj  ,jk))) 
     1531            END_3D 
    15591532         END DO 
    15601533         ! 
    … …  
    15831556      IF( .NOT. sdjf%ln_clim ) THEN    
    15841557                                         WRITE(clname, '(a,"_y",i4.4)' ) TRIM( sdjf%clrootname ), kyear    ! add year 
    1585          IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a, "m",i2.2)' ) TRIM( clname          ), kmonth   ! add month 
     1558         IF( sdjf%clftyp /= 'yearly' )   WRITE(clname, '(a, "m",i2.2)' ) TRIM( clname          ), kmonth   ! add month 
    15861559      ELSE 
    15871560         ! build the new filename if climatological data 
    1588          IF( sdjf%cltype /= 'yearly' )   WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), kmonth   ! add month 
    1589       ENDIF 
    1590       IF(    sdjf%cltype == 'daily' .OR. sdjf%cltype(1:4) == 'week' ) & 
     1561         IF( sdjf%clftyp /= 'yearly' )   WRITE(clname, '(a,"_m",i2.2)' ) TRIM( sdjf%clrootname ), kmonth   ! add month 
     1562      ENDIF 
     1563      IF(    sdjf%clftyp == 'daily' .OR. sdjf%clftyp(1:4) == 'week' ) & 
    15911564         &                               WRITE(clname, '(a,"d" ,i2.2)' ) TRIM( clname          ), kday     ! add day 
    15921565 
    … …  
    16121585         IF( cl_week(ijul) == TRIM(cdday) ) EXIT 
    16131586      END DO 
    1614       IF( ijul .GT. 7 )   CALL ctl_stop( 'ksec_week: wrong day for sdjf%cltype(6:8): '//TRIM(cdday) ) 
     1587      IF( ijul .GT. 7 )   CALL ctl_stop( 'ksec_week: wrong day for sdjf%clftyp(6:8): '//TRIM(cdday) ) 
    16151588      ! 
    16161589      ishift = ijul * NINT(rday) 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/geo2ocean.F90

    r12377 r13540  
    160160      ! (computation done on the north stereographic polar plane) 
    161161      ! 
    162       DO_2D_00_01 
     162      DO_2D( 0, 0, 0, 1 ) 
    163163         !                   
    164164         zlam = plamt(ji,jj)     ! north pole direction & modulous (at t-point) 
    … …  
    249249      ! =============== ! 
    250250 
    251       DO_2D_00_01 
     251      DO_2D( 0, 0, 0, 1 ) 
    252252         IF( MOD( ABS( plamv(ji,jj) - plamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN 
    253253            gsint(ji,jj) = 0. 
    … …  
    272272      ! =========================== ! 
    273273      !           ! lateral boundary cond.: T-, U-, V-, F-pts, sgn 
    274       CALL lbc_lnk_multi( 'geo2ocean', gcost, 'T', -1., gsint, 'T', -1., gcosu, 'U', -1., gsinu, 'U', -1., &  
    275                       &   gcosv, 'V', -1., gsinv, 'V', -1., gcosf, 'F', -1., gsinf, 'F', -1.  ) 
     274      CALL lbc_lnk_multi( 'geo2ocean', gcost, 'T', -1.0_wp, gsint, 'T', -1.0_wp, gcosu, 'U', -1.0_wp, gsinu, 'U', -1.0_wp, &  
     275                      &   gcosv, 'V', -1.0_wp, gsinv, 'V', -1.0_wp, gcosf, 'F', -1.0_wp, gsinf, 'F', -1.0_wp  ) 
    276276      ! 
    277277   END SUBROUTINE angle 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbc_ice.F90

    r12511 r13540  
    6969   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_oce        !: evap - precip over ocean                 [kg/m2/s] 
    7070   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wndm_ice       !: wind speed module at T-point                 [m/s] 
    71    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz         !: wind speed module at T-point                 [m/s] 
     71   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sstfrz         !: sea surface freezing temperature            [degC] 
     72   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rCdU_ice       !: ice-ocean drag at T-point (<0)               [m/s] 
    7273#endif 
    7374 
    … …  
    8990   ! variables used in the coupled interface 
    9091   INTEGER , PUBLIC, PARAMETER ::   jpl = ncat 
    91    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice          ! jpi, jpj 
     92   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice  
    9293    
    9394   ! already defined in ice.F90 for SI3 
    … …  
    9899#endif 
    99100 
    100    REAL(wp), PUBLIC, SAVE ::   cldf_ice = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-] 
     101   REAL(wp), PUBLIC, SAVE ::   pp_cldf = 0.81    !: cloud fraction over sea ice, summer CLIO value   [-] 
    101102 
    102103   !! arrays relating to embedding ice in the ocean 
    … …  
    131132         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce   (jpi,jpj)     ,   & 
    132133         &      qns_oce (jpi,jpj)     , qsr_oce  (jpi,jpj)     , emp_oce    (jpi,jpj)     ,   & 
    133          &      emp_ice (jpi,jpj)     , sstfrz   (jpi,jpj)     , STAT= ierr(2) ) 
     134         &      emp_ice (jpi,jpj)     , sstfrz   (jpi,jpj)     , rCdU_ice   (jpi,jpj)     , STAT= ierr(2) ) 
    134135#endif 
    135136 
    … …  
    167168   LOGICAL         , PUBLIC, PARAMETER ::   lk_si3     = .FALSE.  !: no SI3 ice model 
    168169   LOGICAL         , PUBLIC, PARAMETER ::   lk_cice    = .FALSE.  !: no CICE ice model 
    169    REAL(wp)        , PUBLIC, PARAMETER ::   cldf_ice = 0.81       !: cloud fraction over sea ice, summer CLIO value   [-] 
     170   REAL(wp)        , PUBLIC, PARAMETER ::   pp_cldf    = 0.81     !: cloud fraction over sea ice, summer CLIO value   [-] 
    170171   INTEGER         , PUBLIC, PARAMETER ::   jpl = 1  
    171172   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbc_oce.F90

    r12377 r13540  
    136136   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   atm_co2           !: atmospheric pCO2                             [ppm] 
    137137   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: xcplmask          !: coupling mask for ln_mixcpl (warning: allocated in sbccpl) 
     138   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cloud_fra         !: cloud cover (fraction of cloud in a gridcell) [-] 
    138139 
    139140   !!--------------------------------------------------------------------- 
    … …  
    188189      ! 
    189190      ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     & 
    190          &      atm_co2(jpi,jpj) , tsk_m(jpi,jpj) ,                       & 
     191         &      atm_co2(jpi,jpj) , tsk_m(jpi,jpj) , cloud_fra(jpi,jpj),   & 
    191192         &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) , frq_m(jpi,jpj) ,      & 
    192193         &      ssv_m  (jpi,jpj) , sss_m(jpi,jpj) , ssh_m(jpi,jpj) , STAT=ierr(4) ) 
    … …  
    217218      !!--------------------------------------------------------------------- 
    218219      zcoef = 0.5 / ( zrhoa * zcdrag ) 
    219       DO_2D_00_00 
     220      DO_2D( 0, 0, 0, 0 ) 
    220221         ztx = utau(ji-1,jj  ) + utau(ji,jj) 
    221222         zty = vtau(ji  ,jj-1) + vtau(ji,jj) 
    … …  
    223224         wndm(ji,jj) = SQRT ( ztau * zcoef ) * tmask(ji,jj,1) 
    224225      END_2D 
    225       CALL lbc_lnk( 'sbc_oce', wndm(:,:) , 'T', 1. ) 
     226      CALL lbc_lnk( 'sbc_oce', wndm(:,:) , 'T', 1.0_wp ) 
    226227      ! 
    227228   END SUBROUTINE sbc_tau2wnd 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcapr.F90

    r12511 r13540  
    154154         IF( ln_rstart .AND. iom_varid( numror, 'ssh_ibb', ldstop = .FALSE. ) > 0 ) THEN  
    155155            IF(lwp) WRITE(numout,*) 'sbc_apr:   ssh_ibb read in the restart file' 
    156             CALL iom_get( numror, jpdom_autoglo, 'ssh_ibb', ssh_ibb, ldxios = lrxios )   ! before inv. barometer ssh 
     156            CALL iom_get( numror, jpdom_auto, 'ssh_ibb', ssh_ibb, ldxios = lrxios )   ! before inv. barometer ssh 
    157157            ! 
    158158         ELSE                                         !* no restart: set from nit000 values 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcblk.F90

    r12511 r13540  
    4444   USE lib_fortran    ! to use key_nosignedzero 
    4545#if defined key_si3 
    46    USE ice     , ONLY :   jpl, a_i_b, at_i_b, rn_cnd_s, hfx_err_dif 
    47    USE icethd_dh      ! for CALL ice_thd_snwblow 
     46   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, t_su, rn_cnd_s, hfx_err_dif, nn_qtrice 
     47   USE icevar         ! for CALL ice_var_snwblow 
    4848#endif 
    4949   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009) 
    … …  
    7474#endif 
    7575 
    76    INTEGER , PUBLIC            ::   jpfld         ! maximum number of files to read 
    77    INTEGER , PUBLIC, PARAMETER ::   jp_wndi = 1   ! index of 10m wind velocity (i-component) (m/s)    at T-point 
    78    INTEGER , PUBLIC, PARAMETER ::   jp_wndj = 2   ! index of 10m wind velocity (j-component) (m/s)    at T-point 
    79    INTEGER , PUBLIC, PARAMETER ::   jp_tair = 3   ! index of 10m air temperature             (Kelvin) 
    80    INTEGER , PUBLIC, PARAMETER ::   jp_humi = 4   ! index of specific humidity               ( % ) 
    81    INTEGER , PUBLIC, PARAMETER ::   jp_qsr  = 5   ! index of solar heat                      (W/m2) 
    82    INTEGER , PUBLIC, PARAMETER ::   jp_qlw  = 6   ! index of Long wave                       (W/m2) 
    83    INTEGER , PUBLIC, PARAMETER ::   jp_prec = 7   ! index of total precipitation (rain+snow) (Kg/m2/s) 
    84    INTEGER , PUBLIC, PARAMETER ::   jp_snow = 8   ! index of snow (solid prcipitation)       (kg/m2/s) 
    85    INTEGER , PUBLIC, PARAMETER ::   jp_slp  = 9   ! index of sea level pressure              (Pa) 
    86    INTEGER , PUBLIC, PARAMETER ::   jp_hpgi =10   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
    87    INTEGER , PUBLIC, PARAMETER ::   jp_hpgj =11   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
    88  
     76   INTEGER , PUBLIC, PARAMETER ::   jp_wndi  =  1   ! index of 10m wind velocity (i-component) (m/s)    at T-point 
     77   INTEGER , PUBLIC, PARAMETER ::   jp_wndj  =  2   ! index of 10m wind velocity (j-component) (m/s)    at T-point 
     78   INTEGER , PUBLIC, PARAMETER ::   jp_tair  =  3   ! index of 10m air temperature             (Kelvin) 
     79   INTEGER , PUBLIC, PARAMETER ::   jp_humi  =  4   ! index of specific humidity               ( % ) 
     80   INTEGER , PUBLIC, PARAMETER ::   jp_qsr   =  5   ! index of solar heat                      (W/m2) 
     81   INTEGER , PUBLIC, PARAMETER ::   jp_qlw   =  6   ! index of Long wave                       (W/m2) 
     82   INTEGER , PUBLIC, PARAMETER ::   jp_prec  =  7   ! index of total precipitation (rain+snow) (Kg/m2/s) 
     83   INTEGER , PUBLIC, PARAMETER ::   jp_snow  =  8   ! index of snow (solid prcipitation)       (kg/m2/s) 
     84   INTEGER , PUBLIC, PARAMETER ::   jp_slp   =  9   ! index of sea level pressure              (Pa) 
     85   INTEGER , PUBLIC, PARAMETER ::   jp_uoatm = 10   ! index of surface current (i-component) 
     86   !                                                !          seen by the atmospheric forcing (m/s) at T-point 
     87   INTEGER , PUBLIC, PARAMETER ::   jp_voatm = 11   ! index of surface current (j-component) 
     88   !                                                !          seen by the atmospheric forcing (m/s) at T-point 
     89   INTEGER , PUBLIC, PARAMETER ::   jp_cc    = 12   ! index of cloud cover                     (-)      range:0-1 
     90   INTEGER , PUBLIC, PARAMETER ::   jp_hpgi  = 13   ! index of ABL geostrophic wind or hpg (i-component) (m/s) at T-point 
     91   INTEGER , PUBLIC, PARAMETER ::   jp_hpgj  = 14   ! index of ABL geostrophic wind or hpg (j-component) (m/s) at T-point 
     92   INTEGER , PUBLIC, PARAMETER ::   jpfld    = 14   ! maximum number of files to read 
     93 
     94   ! Warning: keep this structure allocatable for Agrif... 
    8995   TYPE(FLD), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input atmospheric fields (file informations, fields read) 
    9096 
    … …  
    98104   LOGICAL  ::   ln_Cd_L15      ! ice-atm drag = F( ice concentration, atmospheric stability ) (Lupkes et al. JGR2015) 
    99105   ! 
     106   LOGICAL  ::   ln_crt_fbk     ! Add surface current feedback to the wind stress computation  (Renault et al. 2020) 
     107   REAL(wp) ::   rn_stau_a      ! Alpha and Beta coefficients of Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta 
     108   REAL(wp) ::   rn_stau_b      !  
     109   ! 
    100110   REAL(wp)         ::   rn_pfac   ! multiplication factor for precipitation 
    101111   REAL(wp), PUBLIC ::   rn_efac   ! multiplication factor for evaporation 
    102    REAL(wp), PUBLIC ::   rn_vfac   ! multiplication factor for ice/ocean velocity in the calculation of wind stress 
    103112   REAL(wp)         ::   rn_zqt    ! z(q,t) : height of humidity and temperature measurements 
    104113   REAL(wp)         ::   rn_zu     ! z(u)   : height of wind measurements 
    105114   ! 
    106    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice 
    107    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme) 
    108    REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 
     115   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:) ::   Cdn_oce, Chn_oce, Cen_oce  ! neutral coeffs over ocean (L15 bulk scheme and ABL) 
     116   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   Cd_ice , Ch_ice , Ce_ice   ! transfert coefficients over ice 
     117   REAL(wp),         ALLOCATABLE, DIMENSION(:,:) ::   t_zu, q_zu                 ! air temp. and spec. hum. at wind speed height (L15 bulk scheme) 
    109118 
    110119   LOGICAL  ::   ln_skin_cs     ! use the cool-skin (only available in ECMWF and COARE algorithms) !LB 
    … …  
    113122   LOGICAL  ::   ln_humi_dpt    ! humidity read in files ("sn_humi") is dew-point temperature [K] if .true. !LB 
    114123   LOGICAL  ::   ln_humi_rlh    ! humidity read in files ("sn_humi") is relative humidity     [%] if .true. !LB 
     124   LOGICAL  ::   ln_tpot        !!GS: flag to compute or not potential temperature 
    115125   ! 
    116126   INTEGER  ::   nhumi          ! choice of the bulk algorithm 
    … …  
    162172      !! 
    163173      CHARACTER(len=100)            ::   cn_dir                ! Root directory for location of atmospheric forcing files 
    164       TYPE(FLD_N), ALLOCATABLE, DIMENSION(:) ::   slf_i        ! array of namelist informations on the fields to read 
    165       TYPE(FLD_N) ::   sn_wndi, sn_wndj, sn_humi, sn_qsr       ! informations about the fields to be read 
    166       TYPE(FLD_N) ::   sn_qlw , sn_tair, sn_prec, sn_snow      !       "                        " 
    167       TYPE(FLD_N) ::   sn_slp , sn_hpgi, sn_hpgj               !       "                        " 
     174      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                 ! array of namelist informations on the fields to read 
     175      TYPE(FLD_N) ::   sn_wndi, sn_wndj , sn_humi, sn_qsr      ! informations about the fields to be read 
     176      TYPE(FLD_N) ::   sn_qlw , sn_tair , sn_prec, sn_snow     !       "                        " 
     177      TYPE(FLD_N) ::   sn_slp , sn_uoatm, sn_voatm             !       "                        " 
     178      TYPE(FLD_N) ::   sn_cc, sn_hpgi, sn_hpgj                 !       "                        " 
     179      INTEGER     ::   ipka                                    ! number of levels in the atmospheric variable 
    168180      NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields 
    169          &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_hpgi, sn_hpgj,       & 
     181         &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_uoatm, sn_voatm,     & 
     182         &                 sn_cc, sn_hpgi, sn_hpgj,                                   & 
    170183         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p6, ln_ECMWF,             &   ! bulk algorithm 
    171184         &                 cn_dir , rn_zqt, rn_zu,                                    & 
    172          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15,           & 
     185         &                 rn_pfac, rn_efac, ln_Cd_L12, ln_Cd_L15, ln_tpot,           & 
     186         &                 ln_crt_fbk, rn_stau_a, rn_stau_b,                          &   ! current feedback 
    173187         &                 ln_skin_cs, ln_skin_wl, ln_humi_sph, ln_humi_dpt, ln_humi_rlh  ! cool-skin / warm-layer !LB 
    174188      !!--------------------------------------------------------------------- 
    … …  
    242256      !                                   !* set the bulk structure 
    243257      !                                      !- store namelist information in an array 
    244       IF( ln_blk ) jpfld = 9 
    245       IF( ln_abl ) jpfld = 11 
    246       ALLOCATE( slf_i(jpfld) ) 
    247       ! 
    248       slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj 
    249       slf_i(jp_qsr ) = sn_qsr    ;   slf_i(jp_qlw ) = sn_qlw 
    250       slf_i(jp_tair) = sn_tair   ;   slf_i(jp_humi) = sn_humi 
    251       slf_i(jp_prec) = sn_prec   ;   slf_i(jp_snow) = sn_snow 
    252       slf_i(jp_slp ) = sn_slp 
    253       IF( ln_abl ) THEN 
    254          slf_i(jp_hpgi) = sn_hpgi   ;   slf_i(jp_hpgj) = sn_hpgj 
    255       END IF 
     258      ! 
     259      slf_i(jp_wndi ) = sn_wndi    ;   slf_i(jp_wndj ) = sn_wndj 
     260      slf_i(jp_qsr  ) = sn_qsr     ;   slf_i(jp_qlw  ) = sn_qlw 
     261      slf_i(jp_tair ) = sn_tair    ;   slf_i(jp_humi ) = sn_humi 
     262      slf_i(jp_prec ) = sn_prec    ;   slf_i(jp_snow ) = sn_snow 
     263      slf_i(jp_slp  ) = sn_slp     ;   slf_i(jp_cc   ) = sn_cc 
     264      slf_i(jp_uoatm) = sn_uoatm   ;   slf_i(jp_voatm) = sn_voatm 
     265      slf_i(jp_hpgi ) = sn_hpgi    ;   slf_i(jp_hpgj ) = sn_hpgj 
     266      ! 
     267      IF( .NOT. ln_abl ) THEN   ! force to not use jp_hpgi and jp_hpgj, should already be done in namelist_* but we never know... 
     268         slf_i(jp_hpgi)%clname = 'NOT USED' 
     269         slf_i(jp_hpgj)%clname = 'NOT USED' 
     270      ENDIF 
    256271      ! 
    257272      !                                      !- allocate the bulk structure 
    … …  
    264279      DO jfpr= 1, jpfld 
    265280         ! 
    266          IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to zero) 
    267             ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 
    268             sf(jfpr)%fnow(:,:,1) = 0._wp 
     281         IF(   ln_abl    .AND.                                                      & 
     282            &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   & 
     283            &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN 
     284            ipka = jpka   ! ABL: some fields are 3D input 
     285         ELSE 
     286            ipka = 1 
     287         ENDIF 
     288         ! 
     289         ALLOCATE( sf(jfpr)%fnow(jpi,jpj,ipka) ) 
     290         ! 
     291         IF( TRIM(sf(jfpr)%clrootname) == 'NOT USED' ) THEN    !--  not used field  --!   (only now allocated and set to default) 
     292            IF(     jfpr == jp_slp ) THEN 
     293               sf(jfpr)%fnow(:,:,1:ipka) = 101325._wp   ! use standard pressure in Pa 
     294            ELSEIF( jfpr == jp_prec .OR. jfpr == jp_snow .OR. jfpr == jp_uoatm .OR. jfpr == jp_voatm ) THEN 
     295               sf(jfpr)%fnow(:,:,1:ipka) = 0._wp        ! no precip or no snow or no surface currents 
     296            ELSEIF( jfpr == jp_hpgi .OR. jfpr == jp_hpgj ) THEN 
     297               IF( .NOT. ln_abl ) THEN 
     298                  DEALLOCATE( sf(jfpr)%fnow )   ! deallocate as not used in this case 
     299               ELSE 
     300                  sf(jfpr)%fnow(:,:,1:ipka) = 0._wp 
     301               ENDIF 
     302            ELSEIF( jfpr == jp_cc  ) THEN 
     303               sf(jp_cc)%fnow(:,:,1:ipka) = pp_cldf 
     304            ELSE 
     305               WRITE(ctmp1,*) 'sbc_blk_init: no default value defined for field number', jfpr 
     306               CALL ctl_stop( ctmp1 ) 
     307            ENDIF 
    269308         ELSE                                                  !-- used field  --! 
    270             IF(   ln_abl    .AND.                                                      & 
    271                &    ( jfpr == jp_wndi .OR. jfpr == jp_wndj .OR. jfpr == jp_humi .OR.   & 
    272                &      jfpr == jp_hpgi .OR. jfpr == jp_hpgj .OR. jfpr == jp_tair     )  ) THEN   ! ABL: some fields are 3D input 
    273                ALLOCATE( sf(jfpr)%fnow(jpi,jpj,jpka) ) 
    274                IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,jpka,2) ) 
    275             ELSE                                                                                ! others or Bulk fields are 2D fiels 
    276                ALLOCATE( sf(jfpr)%fnow(jpi,jpj,1) ) 
    277                IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,1,2) ) 
    278             ENDIF 
     309            IF( sf(jfpr)%ln_tint )   ALLOCATE( sf(jfpr)%fdta(jpi,jpj,ipka,2) )   ! allocate array for temporal interpolation 
    279310            ! 
    280311            IF( sf(jfpr)%freqh > 0. .AND. MOD( NINT(3600. * sf(jfpr)%freqh), nn_fsbc * NINT(rn_Dt) ) /= 0 )   & 
    281                &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   & 
    282                &                 '               This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' ) 
     312         &  CALL ctl_warn( 'sbc_blk_init: sbcmod timestep rn_Dt*nn_fsbc is NOT a submultiple of atmospheric forcing frequency.',   & 
     313         &                 '               This is not ideal. You should consider changing either rn_Dt or nn_fsbc value...' ) 
    283314         ENDIF 
    284315      END DO 
    … …  
    327358         WRITE(numout,*) '      factor applied on precipitation (total & snow)      rn_pfac      = ', rn_pfac 
    328359         WRITE(numout,*) '      factor applied on evaporation                       rn_efac      = ', rn_efac 
    329          WRITE(numout,*) '      factor applied on ocean/ice velocity                rn_vfac      = ', rn_vfac 
    330360         WRITE(numout,*) '         (form absolute (=0) to relative winds(=1))' 
    331361         WRITE(numout,*) '      use ice-atm drag from Lupkes2012                    ln_Cd_L12    = ', ln_Cd_L12 
    332362         WRITE(numout,*) '      use ice-atm drag from Lupkes2015                    ln_Cd_L15    = ', ln_Cd_L15 
     363         WRITE(numout,*) '      use surface current feedback on wind stress         ln_crt_fbk   = ', ln_crt_fbk 
     364         IF(ln_crt_fbk) THEN 
     365         WRITE(numout,*) '         Renault et al. 2020, eq. 10: Stau = Alpha * Wnd + Beta' 
     366         WRITE(numout,*) '            Alpha                                         rn_stau_a    = ', rn_stau_a 
     367         WRITE(numout,*) '            Beta                                          rn_stau_b    = ', rn_stau_b 
     368         ENDIF 
    333369         ! 
    334370         WRITE(numout,*) 
    … …  
    429465      !                                            ! compute the surface ocean fluxes using bulk formulea 
    430466      IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN 
    431          CALL blk_oce_1( kt, sf(jp_wndi)%fnow(:,:,1), sf(jp_wndj)%fnow(:,:,1),   &   !   <<= in 
    432             &                sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1),   &   !   <<= in 
    433             &                sf(jp_slp )%fnow(:,:,1), sst_m, ssu_m, ssv_m,       &   !   <<= in 
    434             &                sf(jp_qsr )%fnow(:,:,1), sf(jp_qlw )%fnow(:,:,1),   &   !   <<= in (wl/cs) 
    435             &                tsk_m, zssq, zcd_du, zsen, zevp )                       !   =>> out 
    436  
    437          CALL blk_oce_2(     sf(jp_tair)%fnow(:,:,1), sf(jp_qsr )%fnow(:,:,1),   &   !   <<= in 
    438             &                sf(jp_qlw )%fnow(:,:,1), sf(jp_prec)%fnow(:,:,1),   &   !   <<= in 
    439             &                sf(jp_snow)%fnow(:,:,1), tsk_m,                     &   !   <<= in 
    440             &                zsen, zevp )                                            !   <=> in out 
     467         CALL blk_oce_1( kt, sf(jp_wndi )%fnow(:,:,1), sf(jp_wndj )%fnow(:,:,1),   &   !   <<= in 
     468            &                sf(jp_tair )%fnow(:,:,1), sf(jp_humi )%fnow(:,:,1),   &   !   <<= in 
     469            &                sf(jp_slp  )%fnow(:,:,1), sst_m, ssu_m, ssv_m,        &   !   <<= in 
     470            &                sf(jp_uoatm)%fnow(:,:,1), sf(jp_voatm)%fnow(:,:,1),   &   !   <<= in 
     471            &                sf(jp_qsr  )%fnow(:,:,1), sf(jp_qlw  )%fnow(:,:,1),   &   !   <<= in (wl/cs) 
     472            &                tsk_m, zssq, zcd_du, zsen, zevp )                         !   =>> out 
     473 
     474         CALL blk_oce_2(     sf(jp_tair )%fnow(:,:,1), sf(jp_qsr  )%fnow(:,:,1),   &   !   <<= in 
     475            &                sf(jp_qlw  )%fnow(:,:,1), sf(jp_prec )%fnow(:,:,1),   &   !   <<= in 
     476            &                sf(jp_snow )%fnow(:,:,1), tsk_m,                      &   !   <<= in 
     477            &                zsen, zevp )                                              !   <=> in out 
    441478      ENDIF 
    442479      ! 
    … …  
    470507 
    471508 
    472    SUBROUTINE blk_oce_1( kt, pwndi, pwndj , ptair, phumi, &  ! inp 
    473       &                  pslp , pst   , pu   , pv,        &  ! inp 
    474       &                  pqsr , pqlw  ,                   &  ! inp 
    475       &                  ptsk, pssq , pcd_du, psen , pevp   )  ! out 
     509   SUBROUTINE blk_oce_1( kt, pwndi, pwndj, ptair, phumi,        &  ! inp 
     510      &                      pslp , pst  , pu   , pv,            &  ! inp 
     511      &                      puatm, pvatm, pqsr , pqlw ,         &  ! inp 
     512      &                      ptsk , pssq , pcd_du, psen, pevp   )   ! out 
    476513      !!--------------------------------------------------------------------- 
    477514      !!                     ***  ROUTINE blk_oce_1  *** 
    … …  
    498535      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pu     ! surface current at U-point (i-component) [m/s] 
    499536      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pv     ! surface current at V-point (j-component) [m/s] 
     537      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   puatm  ! surface current seen by the atm at T-point (i-component) [m/s] 
     538      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pvatm  ! surface current seen by the atm at T-point (j-component) [m/s] 
    500539      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqsr   ! 
    501540      REAL(wp), INTENT(in   ), DIMENSION(:,:) ::   pqlw   ! 
    … …  
    508547      INTEGER  ::   ji, jj               ! dummy loop indices 
    509548      REAL(wp) ::   zztmp                ! local variable 
     549      REAL(wp) ::   zstmax, zstau 
     550#if defined key_cyclone 
    510551      REAL(wp), DIMENSION(jpi,jpj) ::   zwnd_i, zwnd_j    ! wind speed components at T-point 
     552#endif 
     553      REAL(wp), DIMENSION(jpi,jpj) ::   ztau_i, ztau_j    ! wind stress components at T-point 
    511554      REAL(wp), DIMENSION(jpi,jpj) ::   zU_zu             ! bulk wind speed at height zu  [m/s] 
    512555      REAL(wp), DIMENSION(jpi,jpj) ::   ztpot             ! potential temperature of air at z=rn_zqt [K] 
    … …  
    523566      ptsk(:,:) = pst(:,:) + rt0  ! by default: skin temperature = "bulk SST" (will remain this way if NCAR algorithm used!) 
    524567 
     568      ! --- cloud cover --- ! 
     569      cloud_fra(:,:) = sf(jp_cc)%fnow(:,:,1) 
     570 
    525571      ! ----------------------------------------------------------------------------- ! 
    526572      !      0   Wind components and module at T-point relative to the moving ocean   ! 
    … …  
    532578      zwnd_j(:,:) = 0._wp 
    533579      CALL wnd_cyc( kt, zwnd_i, zwnd_j )    ! add analytical tropical cyclone (Vincent et al. JGR 2012) 
    534       DO_2D_00_00 
    535          pwndi(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
    536          pwndj(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
     580      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     581         zwnd_i(ji,jj) = pwndi(ji,jj) + zwnd_i(ji,jj) 
     582         zwnd_j(ji,jj) = pwndj(ji,jj) + zwnd_j(ji,jj) 
     583         ! ... scalar wind at T-point (not masked) 
     584         wndm(ji,jj) = SQRT( zwnd_i(ji,jj) * zwnd_i(ji,jj) + zwnd_j(ji,jj) * zwnd_j(ji,jj) ) 
     585      END_2D 
     586#else 
     587      ! ... scalar wind module at T-point (not masked) 
     588      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     589         wndm(ji,jj) = SQRT(  pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj)  ) 
    537590      END_2D 
    538591#endif 
    539       DO_2D_00_00 
    540          zwnd_i(ji,jj) = (  pwndi(ji,jj) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  ) 
    541          zwnd_j(ji,jj) = (  pwndj(ji,jj) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  ) 
    542       END_2D 
    543       CALL lbc_lnk_multi( 'sbcblk', zwnd_i, 'T', -1., zwnd_j, 'T', -1. ) 
    544       ! ... scalar wind ( = | U10m - U_oce | ) at T-point (masked) 
    545       wndm(:,:) = SQRT(  zwnd_i(:,:) * zwnd_i(:,:)   & 
    546          &             + zwnd_j(:,:) * zwnd_j(:,:)  ) * tmask(:,:,1) 
    547  
    548592      ! ----------------------------------------------------------------------------- ! 
    549593      !      I   Solar FLUX                                                           ! 
    … …  
    593637         !#LB: because AGRIF hates functions that return something else than a scalar, need to 
    594638         !     use scalar version of gamma_moist() ... 
    595          DO_2D_11_11 
    596             ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 
    597          END_2D 
    598       ENDIF 
    599  
    600  
     639         IF( ln_tpot ) THEN 
     640            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     641               ztpot(ji,jj) = ptair(ji,jj) + gamma_moist( ptair(ji,jj), zqair(ji,jj) ) * rn_zqt 
     642            END_2D 
     643         ELSE 
     644            ztpot = ptair(:,:) 
     645         ENDIF 
     646      ENDIF 
    601647 
    602648      !! Time to call the user-selected bulk parameterization for 
    … …  
    627673 
    628674      END SELECT 
    629  
     675       
     676      IF( iom_use('Cd_oce') )   CALL iom_put("Cd_oce",   zcd_oce * tmask(:,:,1)) 
     677      IF( iom_use('Ce_oce') )   CALL iom_put("Ce_oce",   zce_oce * tmask(:,:,1)) 
     678      IF( iom_use('Ch_oce') )   CALL iom_put("Ch_oce",   zch_oce * tmask(:,:,1)) 
     679      !! LB: mainly here for debugging purpose: 
     680      IF( iom_use('theta_zt') ) CALL iom_put("theta_zt", (ztpot-rt0) * tmask(:,:,1)) ! potential temperature at z=zt 
     681      IF( iom_use('q_zt') )     CALL iom_put("q_zt",     zqair       * tmask(:,:,1)) ! specific humidity       " 
     682      IF( iom_use('theta_zu') ) CALL iom_put("theta_zu", (t_zu -rt0) * tmask(:,:,1)) ! potential temperature at z=zu 
     683      IF( iom_use('q_zu') )     CALL iom_put("q_zu",     q_zu        * tmask(:,:,1)) ! specific humidity       " 
     684      IF( iom_use('ssq') )      CALL iom_put("ssq",      pssq        * tmask(:,:,1)) ! saturation specific humidity at z=0 
     685      IF( iom_use('wspd_blk') ) CALL iom_put("wspd_blk", zU_zu       * tmask(:,:,1)) ! bulk wind speed at z=zu 
     686       
    630687      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    631688         !! ptsk and pssq have been updated!!! 
    … …  
    639696      END IF 
    640697 
    641       !!      CALL iom_put( "Cd_oce", zcd_oce)  ! output value of pure ocean-atm. transfer coef. 
    642       !!      CALL iom_put( "Ch_oce", zch_oce)  ! output value of pure ocean-atm. transfer coef. 
    643  
    644       IF( ABS(rn_zu - rn_zqt) < 0.1_wp ) THEN 
    645          !! If zu == zt, then ensuring once for all that: 
    646          t_zu(:,:) = ztpot(:,:) 
    647          q_zu(:,:) = zqair(:,:) 
    648       ENDIF 
    649  
    650  
    651698      !  Turbulent fluxes over ocean  => BULK_FORMULA @ sbcblk_phy.F90 
    652699      ! ------------------------------------------------------------- 
    653700 
    654701      IF( ln_abl ) THEN         !==  ABL formulation  ==!   multiplication by rho_air and turbulent fluxes computation done in ablstp 
    655          !! FL do we need this multiplication by tmask ... ??? 
    656          DO_2D_11_11 
    657             zztmp = zU_zu(ji,jj) !* tmask(ji,jj,1) 
     702         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     703            zztmp = zU_zu(ji,jj) 
    658704            wndm(ji,jj)   = zztmp                   ! Store zU_zu in wndm to compute ustar2 in ablmod 
    659705            pcd_du(ji,jj) = zztmp * zcd_oce(ji,jj) 
    660706            psen(ji,jj)   = zztmp * zch_oce(ji,jj) 
    661707            pevp(ji,jj)   = zztmp * zce_oce(ji,jj) 
     708            rhoa(ji,jj)   = rho_air( ptair(ji,jj), phumi(ji,jj), pslp(ji,jj) ) 
    662709         END_2D 
    663710      ELSE                      !==  BLK formulation  ==!   turbulent fluxes computation 
    664711         CALL BULK_FORMULA( rn_zu, ptsk(:,:), pssq(:,:), t_zu(:,:), q_zu(:,:), & 
    665             &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),         & 
    666             &               wndm(:,:), zU_zu(:,:), pslp(:,:),                 & 
    667             &               taum(:,:), psen(:,:), zqla(:,:),                  & 
    668             &               pEvap=pevp(:,:), prhoa=rhoa(:,:) ) 
     712            &               zcd_oce(:,:), zch_oce(:,:), zce_oce(:,:),          & 
     713            &               wndm(:,:), zU_zu(:,:), pslp(:,:),                  & 
     714            &               taum(:,:), psen(:,:), zqla(:,:),                   & 
     715            &               pEvap=pevp(:,:), prhoa=rhoa(:,:), pfact_evap=rn_efac ) 
    669716 
    670717         zqla(:,:) = zqla(:,:) * tmask(:,:,1) 
    … …  
    673720         pevp(:,:) = pevp(:,:) * tmask(:,:,1) 
    674721 
    675          ! Tau i and j component on T-grid points, using array "zcd_oce" as a temporary array... 
    676          zcd_oce = 0._wp 
    677          WHERE ( wndm > 0._wp ) zcd_oce = taum / wndm 
    678          zwnd_i = zcd_oce * zwnd_i 
    679          zwnd_j = zcd_oce * zwnd_j 
    680  
    681          CALL iom_put( "taum_oce", taum )   ! output wind stress module 
     722         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     723            IF( wndm(ji,jj) > 0._wp ) THEN 
     724               zztmp = taum(ji,jj) / wndm(ji,jj) 
     725#if defined key_cyclone 
     726               ztau_i(ji,jj) = zztmp * zwnd_i(ji,jj) 
     727               ztau_j(ji,jj) = zztmp * zwnd_j(ji,jj) 
     728#else 
     729               ztau_i(ji,jj) = zztmp * pwndi(ji,jj) 
     730               ztau_j(ji,jj) = zztmp * pwndj(ji,jj) 
     731#endif 
     732            ELSE 
     733               ztau_i(ji,jj) = 0._wp 
     734               ztau_j(ji,jj) = 0._wp                  
     735            ENDIF 
     736         END_2D 
     737 
     738         IF( ln_crt_fbk ) THEN   ! aply eq. 10 and 11 of Renault et al. 2020 (doi: 10.1029/2019MS001715) 
     739            zstmax = MIN( rn_stau_a * 3._wp + rn_stau_b, 0._wp )   ! set the max value of Stau corresponding to a wind of 3 m/s (<0) 
     740            DO_2D( 0, 1, 0, 1 )   ! end at jpj and jpi, as ztau_j(ji,jj+1) ztau_i(ji+1,jj) used in the next loop 
     741               zstau = MIN( rn_stau_a * wndm(ji,jj) + rn_stau_b, zstmax )   ! stau (<0) must be smaller than zstmax 
     742               ztau_i(ji,jj) = ztau_i(ji,jj) + zstau * ( 0.5_wp * ( pu(ji-1,jj  ) + pu(ji,jj) ) - puatm(ji,jj) ) 
     743               ztau_j(ji,jj) = ztau_j(ji,jj) + zstau * ( 0.5_wp * ( pv(ji  ,jj-1) + pv(ji,jj) ) - pvatm(ji,jj) ) 
     744               taum(ji,jj) = SQRT( ztau_i(ji,jj) * ztau_i(ji,jj) + ztau_j(ji,jj) * ztau_j(ji,jj) ) 
     745            END_2D 
     746         ENDIF 
    682747 
    683748         ! ... utau, vtau at U- and V_points, resp. 
    684749         !     Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
    685          !     Note the use of MAX(tmask(i,j),tmask(i+1,j) is to mask tau over ice shelves 
    686          DO_2D_10_10 
    687             utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( zwnd_i(ji,jj) + zwnd_i(ji+1,jj  ) ) & 
    688                &          * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
    689             vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( zwnd_j(ji,jj) + zwnd_j(ji  ,jj+1) ) & 
    690                &          * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
     750         !     Note that coastal wind stress is not used in the code... so this extra care has no effect 
     751         DO_2D( 0, 0, 0, 0 )              ! start loop at 2, in case ln_crt_fbk = T 
     752            utau(ji,jj) = 0.5 * ( 2. - umask(ji,jj,1) ) * ( ztau_i(ji,jj) + ztau_i(ji+1,jj  ) ) & 
     753               &              * MAX(tmask(ji,jj,1),tmask(ji+1,jj,1)) 
     754            vtau(ji,jj) = 0.5 * ( 2. - vmask(ji,jj,1) ) * ( ztau_j(ji,jj) + ztau_j(ji  ,jj+1) ) & 
     755               &              * MAX(tmask(ji,jj,1),tmask(ji,jj+1,1)) 
    691756         END_2D 
    692          CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     757 
     758         IF( ln_crt_fbk ) THEN 
     759            CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1., taum, 'T', -1. ) 
     760         ELSE 
     761            CALL lbc_lnk_multi( 'sbcblk', utau, 'U', -1., vtau, 'V', -1. ) 
     762         ENDIF 
     763 
     764         CALL iom_put( "taum_oce", taum )   ! output wind stress module 
    693765 
    694766         IF(sn_cfctl%l_prtctl) THEN 
    … …  
    766838 
    767839      ! use scalar version of L_vap() for AGRIF compatibility 
    768       DO_2D_11_11 
     840      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    769841         zqla(ji,jj) = - L_vap( ztskk(ji,jj) ) * pevp(ji,jj)    ! Latent Heat flux !!GS: possibility to add a global qla to avoid recomputation after abl update 
    770842      END_2D 
    … …  
    861933      ! 
    862934      INTEGER  ::   ji, jj    ! dummy loop indices 
    863       REAL(wp) ::   zwndi_t , zwndj_t             ! relative wind components at T-point 
    864935      REAL(wp) ::   zootm_su                      ! sea-ice surface mean temperature 
    865936      REAL(wp) ::   zztmp1, zztmp2                ! temporary arrays 
    … …  
    872943      ! ------------------------------------------------------------ ! 
    873944      ! C-grid ice dynamics :   U & V-points (same as ocean) 
    874       DO_2D_00_00 
    875          zwndi_t = (  pwndi(ji,jj) - rn_vfac * 0.5_wp * ( puice(ji-1,jj  ) + puice(ji,jj) )  ) 
    876          zwndj_t = (  pwndj(ji,jj) - rn_vfac * 0.5_wp * ( pvice(ji  ,jj-1) + pvice(ji,jj) )  ) 
    877          wndm_ice(ji,jj) = SQRT( zwndi_t * zwndi_t + zwndj_t * zwndj_t ) * tmask(ji,jj,1) 
     945      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     946         wndm_ice(ji,jj) = SQRT( pwndi(ji,jj) * pwndi(ji,jj) + pwndj(ji,jj) * pwndj(ji,jj) ) 
    878947      END_2D 
    879       CALL lbc_lnk( 'sbcblk', wndm_ice, 'T',  1. ) 
    880948      ! 
    881949      ! Make ice-atm. drag dependent on ice concentration 
    … …  
    888956         Ce_ice(:,:) = Ch_ice(:,:)       ! sensible and latent heat transfer coef. are considered identical 
    889957      ENDIF 
    890  
    891       !! IF ( iom_use("Cd_ice") ) CALL iom_put("Cd_ice", Cd_ice)   ! output value of pure ice-atm. transfer coef. 
    892       !! IF ( iom_use("Ch_ice") ) CALL iom_put("Ch_ice", Ch_ice)   ! output value of pure ice-atm. transfer coef. 
    893  
     958       
     959      IF( iom_use('Cd_ice') ) CALL iom_put("Cd_ice", Cd_ice) 
     960      IF( iom_use('Ce_ice') ) CALL iom_put("Ce_ice", Ce_ice) 
     961      IF( iom_use('Ch_ice') ) CALL iom_put("Ch_ice", Ch_ice) 
     962       
    894963      ! local scalars ( place there for vector optimisation purposes) 
    895       !IF (ln_abl) rhoa  (:,:)  = rho_air( ptair(:,:), phumi(:,:), pslp(:,:) ) !!GS: rhoa must be (re)computed here with ABL to avoid division by zero after (TBI) 
    896964      zcd_dui(:,:) = wndm_ice(:,:) * Cd_ice(:,:) 
    897965 
    898966      IF( ln_blk ) THEN 
    899          ! ------------------------------------------------------------ ! 
    900          !    Wind stress relative to the moving ice ( U10m - U_ice )   ! 
    901          ! ------------------------------------------------------------ ! 
    902          ! C-grid ice dynamics :   U & V-points (same as ocean) 
    903          DO_2D_00_00 
    904             putaui(ji,jj) = 0.5_wp * (  rhoa(ji+1,jj) * zcd_dui(ji+1,jj)             & 
    905                &                      + rhoa(ji  ,jj) * zcd_dui(ji  ,jj)  )          & 
    906                &         * ( 0.5_wp * ( pwndi(ji+1,jj) + pwndi(ji,jj) ) - rn_vfac * puice(ji,jj) ) 
    907             pvtaui(ji,jj) = 0.5_wp * (  rhoa(ji,jj+1) * zcd_dui(ji,jj+1)             & 
    908                &                      + rhoa(ji,jj  ) * zcd_dui(ji,jj  )  )          & 
    909                &         * ( 0.5_wp * ( pwndj(ji,jj+1) + pwndj(ji,jj) ) - rn_vfac * pvice(ji,jj) ) 
     967         ! ---------------------------------------------------- ! 
     968         !    Wind stress relative to nonmoving ice ( U10m )    ! 
     969         ! ---------------------------------------------------- ! 
     970         ! supress moving ice in wind stress computation as we don't know how to do it properly... 
     971         DO_2D( 0, 1, 0, 1 )    ! at T point  
     972            putaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndi(ji,jj) 
     973            pvtaui(ji,jj) = rhoa(ji,jj) * zcd_dui(ji,jj) * pwndj(ji,jj) 
    910974         END_2D 
    911          CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1., pvtaui, 'V', -1. ) 
     975         ! 
     976         DO_2D( 0, 0, 0, 0 )    ! U & V-points (same as ocean). 
     977            ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology  
     978            zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) ) 
     979            zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) ) 
     980            putaui(ji,jj) = zztmp1 * ( putaui(ji,jj) + putaui(ji+1,jj  ) ) 
     981            pvtaui(ji,jj) = zztmp2 * ( pvtaui(ji,jj) + pvtaui(ji  ,jj+1) ) 
     982         END_2D 
     983         CALL lbc_lnk_multi( 'sbcblk', putaui, 'U', -1._wp, pvtaui, 'V', -1._wp ) 
    912984         ! 
    913985         IF(sn_cfctl%l_prtctl)  CALL prt_ctl( tab2d_1=putaui  , clinfo1=' blk_ice: putaui : '   & 
    914986            &                               , tab2d_2=pvtaui  , clinfo2='          pvtaui : ' ) 
    915       ELSE 
     987      ELSE ! ln_abl 
    916988         zztmp1 = 11637800.0_wp 
    917989         zztmp2 =    -5897.8_wp 
    918          DO_2D_11_11 
     990         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    919991            pcd_dui(ji,jj) = zcd_dui (ji,jj) 
    920992            pseni  (ji,jj) = wndm_ice(ji,jj) * Ch_ice(ji,jj) 
    … …  
    9571029      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    9581030      REAL(wp) ::   zztmp, zztmp2, z1_rLsub  !   -      - 
    959       REAL(wp) ::   zfr1, zfr2               ! local variables 
    9601031      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_st         ! inverse of surface temperature 
    9611032      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_qlw         ! long wave heat flux over ice 
    … …  
    9661037      REAL(wp), DIMENSION(jpi,jpj)     ::   zqair         ! specific humidity of air at z=rn_zqt [kg/kg] !LB 
    9671038      REAL(wp), DIMENSION(jpi,jpj)     ::   ztmp, ztmp2 
     1039      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri 
    9681040      !!--------------------------------------------------------------------- 
    9691041      ! 
    … …  
    10461118      evap_ice (:,:,:) = rn_efac * qla_ice (:,:,:) * z1_rLsub    ! sublimation 
    10471119      devap_ice(:,:,:) = rn_efac * dqla_ice(:,:,:) * z1_rLsub    ! d(sublimation)/dT 
    1048       zevap    (:,:)   = rn_efac * ( emp(:,:) + tprecip(:,:) )   ! evaporation over ocean 
     1120      zevap    (:,:)   = emp(:,:) + tprecip(:,:)   ! evaporation over ocean  !LB: removed rn_efac here, correct??? 
    10491121 
    10501122      ! --- evaporation minus precipitation --- ! 
    10511123      zsnw(:,:) = 0._wp 
    1052       CALL ice_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing 
     1124      CALL ice_var_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing 
    10531125      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw ) 
    10541126      emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw 
    … …  
    10771149      END DO 
    10781150 
    1079       ! --- shortwave radiation transmitted below the surface (W/m2, see Grenfell Maykut 77) --- ! 
    1080       zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            ! transmission when hi>10cm 
    1081       zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            ! zfr2 such that zfr1 + zfr2 to equal 1 
    1082       ! 
    1083       WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm 
    1084          qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( zfr1 + zfr2 * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
    1085       ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (zfr1) when hi>10cm 
    1086          qtr_ice_top(:,:,:) = qsr_ice(:,:,:) * zfr1 
    1087       ELSEWHERE                                                         ! zero when hs>0 
    1088          qtr_ice_top(:,:,:) = 0._wp 
    1089       END WHERE 
    1090       ! 
    1091  
     1151      ! --- shortwave radiation transmitted thru the surface scattering layer (W/m2) --- ! 
     1152      IF( nn_qtrice == 0 ) THEN 
     1153         ! formulation derived from Grenfell and Maykut (1977), where transmission rate 
     1154         !    1) depends on cloudiness 
     1155         !    2) is 0 when there is any snow 
     1156         !    3) tends to 1 for thin ice 
     1157         ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm 
     1158         DO jl = 1, jpl 
     1159            WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )     ! linear decrease from hi=0 to 10cm   
     1160               qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) ) 
     1161            ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )     ! constant (ztri) when hi>10cm 
     1162               qtr_ice_top(:,:,jl) = qsr_ice(:,:,jl) * ztri(:,:) 
     1163            ELSEWHERE                                                         ! zero when hs>0 
     1164               qtr_ice_top(:,:,jl) = 0._wp  
     1165            END WHERE 
     1166         ENDDO 
     1167      ELSEIF( nn_qtrice == 1 ) THEN 
     1168         ! formulation is derived from the thesis of M. Lebrun (2019). 
     1169         !    It represents the best fit using several sets of observations 
     1170         !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90) 
     1171         qtr_ice_top(:,:,:) = 0.3_wp * qsr_ice(:,:,:) 
     1172      ENDIF 
     1173      ! 
    10921174      IF( iom_use('evap_ao_cea') .OR. iom_use('hflx_evap_cea') ) THEN 
    10931175         ztmp(:,:) = zevap(:,:) * ( 1._wp - at_i_b(:,:) ) 
    … …  
    11711253         ! 
    11721254         DO jl = 1, jpl 
    1173             DO_2D_11_11 
     1255            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    11741256               zhe = ( rn_cnd_s * phi(ji,jj,jl) + rcnd_i * phs(ji,jj,jl) ) * zfac                            ! Effective thickness 
    11751257               IF( zhe >=  zfac2 )   zgfac(ji,jj,jl) = MIN( 2._wp, 0.5_wp * ( 1._wp + LOG( zhe * zfac3 ) ) ) ! Enhanced conduction factor 
    … …  
    11861268      ! 
    11871269      DO jl = 1, jpl 
    1188          DO_2D_11_11 
     1270         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    11891271            ! 
    11901272            zkeff_h = zfac * zgfac(ji,jj,jl) / &                                    ! Effective conductivity of the snow-ice system divided by thickness 
    … …  
    13341416      zqi_sat(:,:) =                  q_sat( ptm_su(:,:), pslp(:,:) )   ! saturation humidity over ice   [kg/kg] 
    13351417      ! 
    1336       DO_2D_00_00 
     1418      DO_2D( 0, 0, 0, 0 ) 
    13371419         ! Virtual potential temperature [K] 
    13381420         zthetav_os = zst(ji,jj)    * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean 
    … …  
    13771459         ! 
    13781460      END_2D 
    1379       CALL lbc_lnk_multi( 'sbcblk', pcd, 'T',  1., pch, 'T', 1. ) 
     1461      CALL lbc_lnk_multi( 'sbcblk', pcd, 'T',  1.0_wp, pch, 'T', 1.0_wp ) 
    13801462      ! 
    13811463   END SUBROUTINE Cdn10_Lupkes2015 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcblk_algo_coare3p0.F90

    r12377 r13540  
    194194      IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P0_INIT(l_use_cs, l_use_wl) 
    195195 
    196       l_zt_equal_zu = .FALSE. 
    197       IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     196      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 
    198197      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) ) 
    199198 
    … …  
    395394      !!------------------------------------------------------------------- 
    396395      ! 
    397       DO_2D_11_11 
    398          ! 
    399          zw = pwnd(ji,jj)   ! wind speed 
    400          ! 
    401          ! Charnock's constant, increases with the wind : 
    402          zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10))  ! If zw<10. --> 0, else --> 1 
    403          zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1 
    404          ! 
    405          alfa_charn_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s 
    406             &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) & 
    407             &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999) 
    408          ! 
     396      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     397      ! 
     398      zw = pwnd(ji,jj)   ! wind speed 
     399      ! 
     400      ! Charnock's constant, increases with the wind : 
     401      zgt10 = 0.5 + SIGN(0.5_wp,(zw - 10))  ! If zw<10. --> 0, else --> 1 
     402      zgt18 = 0.5 + SIGN(0.5_wp,(zw - 18.)) ! If zw<18. --> 0, else --> 1 
     403      ! 
     404      alfa_charn_3p0(ji,jj) =  (1. - zgt10)*0.011    &    ! wind is lower than 10 m/s 
     405         &     + zgt10*((1. - zgt18)*(0.011 + (0.018 - 0.011) & 
     406         &      *(zw - 10.)/(18. - 10.)) + zgt18*( 0.018 ) )    ! Hare et al. (1999) 
     407      ! 
    409408      END_2D 
    410409      ! 
    … …  
    431430      !!---------------------------------------------------------------------------------- 
    432431      ! 
    433       DO_2D_11_11 
    434          ! 
    435          zta = pzeta(ji,jj) 
    436          ! 
    437          zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable 
    438          ! 
    439          zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   & 
    440             & - 2.*ATAN(zphi_m) + 0.5*rpi 
    441          ! 
    442          zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective 
    443          ! 
    444          zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
    445             &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
    446          ! 
    447          zf = zta*zta 
    448          zf = zf/(1. + zf) 
    449          zc = MIN(50._wp, 0.35_wp*zta) 
    450          zstab = 0.5 + SIGN(0.5_wp, zta) 
    451          ! 
    452          psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 
    453             &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0) 
    454             &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     " 
    455          ! 
     432      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     433      ! 
     434      zta = pzeta(ji,jj) 
     435      ! 
     436      zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable 
     437      ! 
     438      zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   & 
     439         & - 2.*ATAN(zphi_m) + 0.5*rpi 
     440      ! 
     441      zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective 
     442      ! 
     443      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
     444         &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
     445      ! 
     446      zf = zta*zta 
     447      zf = zf/(1. + zf) 
     448      zc = MIN(50._wp, 0.35_wp*zta) 
     449      zstab = 0.5 + SIGN(0.5_wp, zta) 
     450      ! 
     451      psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 
     452         &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0) 
     453         &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     " 
     454      ! 
    456455      END_2D 
    457456      ! 
    … …  
    482481      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab 
    483482      ! 
    484       DO_2D_11_11 
    485          ! 
    486          zta = pzeta(ji,jj) 
    487          ! 
    488          zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 
    489          ! 
    490          zpsi_k = 2.*LOG((1. + zphi_h)/2.) 
    491          ! 
    492          zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective 
    493          ! 
    494          zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
    495             &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
    496          ! 
    497          zf = zta*zta 
    498          zf = zf/(1. + zf) 
    499          zc = MIN(50._wp,0.35_wp*zta) 
    500          zstab = 0.5 + SIGN(0.5_wp, zta) 
    501          ! 
    502          psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 
    503             &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     & 
    504             &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 
    505          ! 
     483      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     484      ! 
     485      zta = pzeta(ji,jj) 
     486      ! 
     487      zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 
     488      ! 
     489      zpsi_k = 2.*LOG((1. + zphi_h)/2.) 
     490      ! 
     491      zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective 
     492      ! 
     493      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
     494         &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
     495      ! 
     496      zf = zta*zta 
     497      zf = zf/(1. + zf) 
     498      zc = MIN(50._wp,0.35_wp*zta) 
     499      zstab = 0.5 + SIGN(0.5_wp, zta) 
     500      ! 
     501      psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 
     502         &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     & 
     503         &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 
     504      ! 
    506505      END_2D 
    507506      ! 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcblk_algo_coare3p6.F90

    r12377 r13540  
    194194      IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl) 
    195195 
    196       l_zt_equal_zu = .FALSE. 
    197       IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     196      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 
    198197      IF( .NOT. l_zt_equal_zu )  ALLOCATE( zeta_t(jpi,jpj) ) 
    199198 
    … …  
    431430      !!---------------------------------------------------------------------------------- 
    432431      ! 
    433       DO_2D_11_11 
    434          ! 
    435          zta = pzeta(ji,jj) 
    436          ! 
    437          zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable 
    438          ! 
    439          zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   & 
    440             & - 2.*ATAN(zphi_m) + 0.5*rpi 
    441          ! 
    442          zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective 
    443          ! 
    444          zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
    445             &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
    446          ! 
    447          zf = zta*zta 
    448          zf = zf/(1. + zf) 
    449          zc = MIN(50._wp, 0.35_wp*zta) 
    450          zstab = 0.5 + SIGN(0.5_wp, zta) 
    451          ! 
    452          psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 
    453             &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0) 
    454             &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     " 
    455          ! 
     432      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     433      ! 
     434      zta = pzeta(ji,jj) 
     435      ! 
     436      zphi_m = ABS(1. - 15.*zta)**.25    !!Kansas unstable 
     437      ! 
     438      zpsi_k = 2.*LOG((1. + zphi_m)/2.) + LOG((1. + zphi_m*zphi_m)/2.)   & 
     439         & - 2.*ATAN(zphi_m) + 0.5*rpi 
     440      ! 
     441      zphi_c = ABS(1. - 10.15*zta)**.3333                   !!Convective 
     442      ! 
     443      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
     444         &     - 1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
     445      ! 
     446      zf = zta*zta 
     447      zf = zf/(1. + zf) 
     448      zc = MIN(50._wp, 0.35_wp*zta) 
     449      zstab = 0.5 + SIGN(0.5_wp, zta) 
     450      ! 
     451      psi_m_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & ! (zta < 0) 
     452         &                -   zstab     * ( 1. + 1.*zta     &                ! (zta > 0) 
     453         &                         + 0.6667*(zta - 14.28)/EXP(zc) + 8.525 )   !     " 
     454      ! 
    456455      END_2D 
    457456      ! 
    … …  
    482481      REAL(wp) :: zta, zphi_h, zphi_c, zpsi_k, zpsi_c, zf, zc, zstab 
    483482      ! 
    484       DO_2D_11_11 
    485          ! 
    486          zta = pzeta(ji,jj) 
    487          ! 
    488          zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 
    489          ! 
    490          zpsi_k = 2.*LOG((1. + zphi_h)/2.) 
    491          ! 
    492          zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective 
    493          ! 
    494          zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
    495             &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
    496          ! 
    497          zf = zta*zta 
    498          zf = zf/(1. + zf) 
    499          zc = MIN(50._wp,0.35_wp*zta) 
    500          zstab = 0.5 + SIGN(0.5_wp, zta) 
    501          ! 
    502          psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 
    503             &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     & 
    504             &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 
    505          ! 
     483      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     484      ! 
     485      zta = pzeta(ji,jj) 
     486      ! 
     487      zphi_h = (ABS(1. - 15.*zta))**.5  !! Kansas unstable   (zphi_h = zphi_m**2 when unstable, zphi_m when stable) 
     488      ! 
     489      zpsi_k = 2.*LOG((1. + zphi_h)/2.) 
     490      ! 
     491      zphi_c = (ABS(1. - 34.15*zta))**.3333   !! Convective 
     492      ! 
     493      zpsi_c = 1.5*LOG((1. + zphi_c + zphi_c*zphi_c)/3.) & 
     494         &    -1.7320508*ATAN((1. + 2.*zphi_c)/1.7320508) + 1.813799447 
     495      ! 
     496      zf = zta*zta 
     497      zf = zf/(1. + zf) 
     498      zc = MIN(50._wp,0.35_wp*zta) 
     499      zstab = 0.5 + SIGN(0.5_wp, zta) 
     500      ! 
     501      psi_h_coare(ji,jj) = (1. - zstab) * ( (1. - zf)*zpsi_k + zf*zpsi_c ) & 
     502         &                -   zstab     * ( (ABS(1. + 2.*zta/3.))**1.5     & 
     503         &                           + .6667*(zta - 14.28)/EXP(zc) + 8.525 ) 
     504      ! 
    506505      END_2D 
    507506      ! 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r12377 r13540  
    9898      &                      Qsw, rad_lw, slp, pdT_cs,                                & ! optionals for cool-skin (and warm-layer) 
    9999      &                      pdT_wl, pHz_wl )                                           ! optionals for warm-layer only 
    100       !!---------------------------------------------------------------------- 
     100      !!---------------------------------------------------------------------------------- 
    101101      !!                      ***  ROUTINE  turb_ecmwf  *** 
    102102      !! 
    … …  
    184184      LOGICAL :: l_zt_equal_zu = .FALSE.      ! if q and t are given at same height as U 
    185185      ! 
    186       REAL(wp), DIMENSION(jpi,jpj) ::  u_star, t_star, q_star 
    187       REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu      
    188       REAL(wp), DIMENSION(jpi,jpj) :: znu_a !: Nu_air, Viscosity of air 
     186      REAL(wp), DIMENSION(jpi,jpj) :: u_star, t_star, q_star 
     187      REAL(wp), DIMENSION(jpi,jpj) :: dt_zu, dq_zu 
     188      REAL(wp), DIMENSION(jpi,jpj) :: znu_a         !: Nu_air, Viscosity of air 
    189189      REAL(wp), DIMENSION(jpi,jpj) :: Linv  !: 1/L (inverse of Monin Obukhov length... 
    190190      REAL(wp), DIMENSION(jpi,jpj) :: z0, z0t, z0q 
    … …  
    196196      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_ecmwf@sbcblk_algo_ecmwf.F90' 
    197197      !!---------------------------------------------------------------------------------- 
    198  
    199198      IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 
    200199 
    201       l_zt_equal_zu = .FALSE. 
    202       IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     200      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 
    203201 
    204202      !! Initializations for cool skin and warm layer: 
    … …  
    412410      REAL(wp) :: zzeta, zx, ztmp, psi_unst, psi_stab, stab 
    413411      !!---------------------------------------------------------------------------------- 
    414       DO_2D_11_11 
    415          ! 
    416          zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 
    417          ! 
    418          ! Unstable (Paulson 1970): 
    419          !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    420          zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 
    421          ztmp = 1._wp + SQRT(zx) 
    422          ztmp = ztmp*ztmp 
    423          psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   & 
    424             &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 
    425          ! 
    426          ! Unstable: 
    427          ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    428          psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 
    429             &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp 
    430          ! 
    431          ! Combining: 
    432          stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
    433          ! 
    434          psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 
    435             &                +      stab  * psi_stab      ! (zzeta > 0) Stable 
    436          ! 
     412      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     413      ! 
     414      zzeta = MIN( pzeta(ji,jj) , 5._wp ) !! Very stable conditions (L positif and big!): 
     415      ! 
     416      ! Unstable (Paulson 1970): 
     417      !   eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
     418      zx = SQRT(ABS(1._wp - 16._wp*zzeta)) 
     419      ztmp = 1._wp + SQRT(zx) 
     420      ztmp = ztmp*ztmp 
     421      psi_unst = LOG( 0.125_wp*ztmp*(1._wp + zx) )   & 
     422         &       -2._wp*ATAN( SQRT(zx) ) + 0.5_wp*rpi 
     423      ! 
     424      ! Unstable: 
     425      ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
     426      psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & 
     427         &       - zzeta - 2._wp/3._wp*5._wp/0.35_wp 
     428      ! 
     429      ! Combining: 
     430      stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
     431      ! 
     432      psi_m_ecmwf(ji,jj) = (1._wp - stab) * psi_unst & ! (zzeta < 0) Unstable 
     433         &                +      stab  * psi_stab      ! (zzeta > 0) Stable 
     434      ! 
    437435      END_2D 
    438436   END FUNCTION psi_m_ecmwf 
    … …  
    457455      !!---------------------------------------------------------------------------------- 
    458456      ! 
    459       DO_2D_11_11 
    460          ! 
    461          zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!): 
    462          ! 
    463          zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!! 
    464          !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 
    465          ! Unstable (Paulson 1970) : 
    466          psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
    467          ! 
    468          ! Stable: 
    469          psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
    470             &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 
    471          ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
    472          ! 
    473          stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
    474          ! 
    475          ! 
    476          psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable 
    477             &                +    stab    * psi_stab        ! (zzeta > 0) Stable 
    478          ! 
     457      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     458      ! 
     459      zzeta = MIN(pzeta(ji,jj) , 5._wp)   ! Very stable conditions (L positif and big!): 
     460      ! 
     461      zx  = ABS(1._wp - 16._wp*zzeta)**.25        ! this is actually (1/phi_m)**2  !!! 
     462      !                                     ! eq.3.19, Chap.3, p.33, IFS doc - Cy31r1 
     463      ! Unstable (Paulson 1970) : 
     464      psi_unst = 2._wp*LOG(0.5_wp*(1._wp + zx*zx))   ! eq.3.20, Chap.3, p.33, IFS doc - Cy31r1 
     465      ! 
     466      ! Stable: 
     467      psi_stab = -2._wp/3._wp*(zzeta - 5._wp/0.35_wp)*EXP(-0.35_wp*zzeta) & ! eq.3.22, Chap.3, p.33, IFS doc - Cy31r1 
     468         &       - ABS(1._wp + 2._wp/3._wp*zzeta)**1.5_wp - 2._wp/3._wp*5._wp/0.35_wp + 1._wp 
     469      ! LB: added ABS() to avoid NaN values when unstable, which contaminates the unstable solution... 
     470      ! 
     471      stab = 0.5_wp + SIGN(0.5_wp, zzeta) ! zzeta > 0 => stab = 1 
     472      ! 
     473      ! 
     474      psi_h_ecmwf(ji,jj) = (1._wp - stab) * psi_unst &   ! (zzeta < 0) Unstable 
     475         &                +    stab    * psi_stab        ! (zzeta > 0) Stable 
     476      ! 
    479477      END_2D 
    480478   END FUNCTION psi_h_ecmwf 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcblk_algo_ncar.F90

    r12377 r13540  
    112112      REAL(wp), DIMENSION(jpi,jpj) ::   stab          ! stability test integer 
    113113      !!---------------------------------------------------------------------------------- 
    114       ! 
    115       l_zt_equal_zu = .FALSE. 
    116       IF( ABS(zu - zt) < 0.01_wp )   l_zt_equal_zu = .TRUE.    ! testing "zu == zt" is risky with double precision 
     114      l_zt_equal_zu = ( ABS(zu - zt) < 0.01_wp ) ! testing "zu == zt" is risky with double precision 
    117115 
    118116      U_blk = MAX( 0.5_wp , U_zu )   !  relative wind speed at zu (normally 10m), we don't want to fall under 0.5 m/s 
    … …  
    143141      ENDIF 
    144142 
    145       !! Initializing values at z_u with z_t values: 
    146       t_zu = t_zt   ;   q_zu = q_zt 
     143      !! First guess of temperature and humidity at height zu: 
     144      t_zu = MAX( t_zt ,  180._wp )   ! who knows what's given on masked-continental regions... 
     145      q_zu = MAX( q_zt , 1.e-6_wp )   !               " 
    147146 
    148147      !! ITERATION BLOCK 
    … …  
    242241      !!---------------------------------------------------------------------------------- 
    243242      ! 
    244       DO_2D_11_11 
     243      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    245244         ! 
    246245         zw  = pw10(ji,jj) 
    … …  
    278277      REAL(wp) :: zx2, zx, zstab   ! local scalars 
    279278      !!---------------------------------------------------------------------------------- 
    280       DO_2D_11_11 
     279      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    281280         zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 
    282281         zx2 = MAX( zx2 , 1._wp ) 
    … …  
    309308      !!---------------------------------------------------------------------------------- 
    310309      ! 
    311       DO_2D_11_11 
     310      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    312311         zx2 = SQRT( ABS( 1._wp - 16._wp*pzeta(ji,jj) ) ) 
    313312         zx2 = MAX( zx2 , 1._wp ) 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcblk_phy.F90

    r12377 r13540  
    3131   REAL(wp), PARAMETER, PUBLIC :: R_vap   = 461.495_wp  !: Specific gas constant for water vapor          [J/K/kg] 
    3232   REAL(wp), PARAMETER, PUBLIC :: reps0   = R_dry/R_vap !: ratio of gas constant for dry air and water vapor => ~ 0.622 
    33    REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
     33   REAL(wp), PARAMETER, PUBLIC :: rctv0   = R_vap/R_dry - 1._wp  !: for virtual temperature (== (1-eps)/eps) => ~ 0.608 
    3434   REAL(wp), PARAMETER, PUBLIC :: rCp_air = 1000.5_wp   !: specific heat of air (only used for ice fluxes now...) 
    3535   REAL(wp), PARAMETER, PUBLIC :: rCd_ice = 1.4e-3_wp   !: transfer coefficient over ice 
    … …  
    181181      !!---------------------------------------------------------------------------------- 
    182182      ! 
    183       DO_2D_11_11 
     183      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    184184         ztc  = ptak(ji,jj) - rt0   ! air temp, in deg. C 
    185185         ztc2 = ztc*ztc 
    … …  
    270270      INTEGER  ::   ji, jj         ! dummy loop indices 
    271271      !!---------------------------------------------------------------------------------- 
    272       DO_2D_11_11 
     272      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    273273         gamma_moist_vctr(ji,jj) = gamma_moist_sclr( ptak(ji,jj), pqa(ji,jj) ) 
    274274      END_2D 
    … …  
    315315      !!------------------------------------------------------------------- 
    316316      ! 
    317       DO_2D_11_11 
     317      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    318318         ! 
    319319         zqa = (1._wp + rctv0*pqa(ji,jj)) 
    … …  
    351351      !!------------------------------------------------------------------- 
    352352      ! 
    353       DO_2D_11_11 
     353      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    354354         ! 
    355355         zqa = 0.5_wp*(pqa(ji,jj)+pssq(ji,jj))                                        ! ~ mean q within the layer... 
    … …  
    448448      !!---------------------------------------------------------------------------------- 
    449449      ! 
    450       DO_2D_11_11 
     450      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    451451         ! 
    452452         ze_sat =  e_sat_sclr( ptak(ji,jj) ) 
    … …  
    473473      !!---------------------------------------------------------------------------------- 
    474474      ! 
    475       DO_2D_11_11 
     475      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    476476         ze = prha(ji,jj)*e_sat_sclr(ptak(ji,jj)) 
    477477         q_air_rh(ji,jj) = ze*reps0/(pslp(ji,jj) - (1. - reps0)*ze) 
    … …  
    511511      INTEGER  ::   ji, jj     ! dummy loop indices 
    512512      !!---------------------------------------------------------------------------------- 
    513       DO_2D_11_11 
     513      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    514514 
    515515         zdt = pTa(ji,jj) - pTs(ji,jj) ;  zdt = SIGN( MAX(ABS(zdt),1.E-6_wp), zdt ) 
    … …  
    520520         zCe = zz0*pqst(ji,jj)/zdq 
    521521 
    522          CALL BULK_FORMULA( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), & 
    523             &              zCd, zCh, zCe,                                        & 
    524             &              pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                 & 
    525             &              pTau(ji,jj), zQsen, zQlat ) 
    526  
     522         CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), & 
     523            &                    zCd, zCh, zCe,                                       & 
     524            &                    pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                & 
     525            &                    pTau(ji,jj), zQsen, zQlat ) 
     526          
    527527         zTs2  = pTs(ji,jj)*pTs(ji,jj) 
    528528         zQlw  = emiss_w*(prlw(ji,jj) - stefan*zTs2*zTs2) ! Net longwave flux 
    … …  
    535535 
    536536 
    537    SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa,  & 
    538       &                          pCd, pCh, pCe,            & 
    539       &                          pwnd, pUb, pslp,          & 
    540       &                          pTau, pQsen, pQlat,  pEvap, prhoa ) 
     537   SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, & 
     538      &                          pCd, pCh, pCe,           & 
     539      &                          pwnd, pUb, pslp,         & 
     540      &                          pTau, pQsen, pQlat,      & 
     541      &                          pEvap, prhoa, pfact_evap ) 
     542      !!---------------------------------------------------------------------------------- 
     543      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m) 
     544      REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K] 
     545      REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg] 
     546      REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K] 
     547      REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg] 
     548      REAL(wp), INTENT(in)  :: pCd 
     549      REAL(wp), INTENT(in)  :: pCh 
     550      REAL(wp), INTENT(in)  :: pCe 
     551      REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s] 
     552      REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 
     553      REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa] 
     554      !! 
     555      REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2] 
     556      REAL(wp), INTENT(out) :: pQsen !  [W/m^2] 
     557      REAL(wp), INTENT(out) :: pQlat !  [W/m^2] 
     558      !! 
     559      REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 
     560      REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3] 
     561      REAL(wp), INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent) 
     562      !! 
     563      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap 
     564      INTEGER  :: jq 
     565      !!---------------------------------------------------------------------------------- 
     566      zfact_evap = 1._wp 
     567      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap 
     568       
     569      !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 
     570      ztaa = pTa ! first guess... 
     571      DO jq = 1, 4 
     572         zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa )  !LOLO: why not "0.5*(pqs+pqa)" rather then "pqa" ??? 
     573         ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder... 
     574      END DO 
     575      zrho = rho_air(ztaa, pqa, pslp) 
     576      zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 
     577 
     578      zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10 
     579 
     580      pTau = zUrho * pCd * pwnd ! Wind stress module 
     581 
     582      zevap = zUrho * pCe * (pqa - pqs) 
     583      pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa) 
     584      pQlat = L_vap(pTs) * zevap 
     585 
     586      IF( PRESENT(pEvap) ) pEvap = - zfact_evap * zevap 
     587      IF( PRESENT(prhoa) ) prhoa = zrho 
     588 
     589   END SUBROUTINE BULK_FORMULA_SCLR 
     590 
     591   SUBROUTINE BULK_FORMULA_VCTR( pzu, pTs, pqs, pTa, pqa, & 
     592      &                          pCd, pCh, pCe,           & 
     593      &                          pwnd, pUb, pslp,         & 
     594      &                          pTau, pQsen, pQlat,      &  
     595      &                          pEvap, prhoa, pfact_evap )       
    541596      !!---------------------------------------------------------------------------------- 
    542597      REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m) 
    … …  
    558613      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 
    559614      REAL(wp), DIMENSION(jpi,jpj), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3] 
    560       !! 
    561       REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap 
    562       INTEGER  :: ji, jj, jq     ! dummy loop indices 
    563       !!---------------------------------------------------------------------------------- 
    564       DO_2D_11_11 
    565  
    566          !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 
    567          ztaa = pTa(ji,jj) ! first guess... 
    568          DO jq = 1, 4 
    569             zgamma = gamma_moist( 0.5*(ztaa+pTs(ji,jj)) , pqa(ji,jj) ) 
    570             ztaa = pTa(ji,jj) - zgamma*pzu   ! Absolute temp. is slightly colder... 
    571          END DO 
    572          zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)) 
    573          zrho = rho_air(ztaa, pqa(ji,jj), pslp(ji,jj)-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 
    574  
    575          zUrho = pUb(ji,jj)*MAX(zrho, 1._wp)     ! rho*U10 
    576  
    577          pTau(ji,jj) = zUrho * pCd(ji,jj) * pwnd(ji,jj) ! Wind stress module 
    578  
    579          zevap        = zUrho * pCe(ji,jj) * (pqa(ji,jj) - pqs(ji,jj)) 
    580          pQsen(ji,jj) = zUrho * pCh(ji,jj) * (pTa(ji,jj) - pTs(ji,jj)) * cp_air(pqa(ji,jj)) 
    581          pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap 
    582  
    583          IF( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 
     615      REAL(wp),                     INTENT(in) , OPTIONAL :: pfact_evap  ! ABOMINATION: corrective factor for evaporation (doing this against my will! /laurent) 
     616      !! 
     617      REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap, zfact_evap 
     618      INTEGER  :: ji, jj 
     619      !!---------------------------------------------------------------------------------- 
     620      zfact_evap = 1._wp 
     621      IF( PRESENT(pfact_evap) ) zfact_evap = pfact_evap 
     622 
     623      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     624 
     625         CALL BULK_FORMULA_SCLR( pzu, pTs(ji,jj), pqs(ji,jj), pTa(ji,jj), pqa(ji,jj), & 
     626            &                    pCd(ji,jj), pCh(ji,jj), pCe(ji,jj),                  & 
     627            &                    pwnd(ji,jj), pUb(ji,jj), pslp(ji,jj),                & 
     628            &                    pTau(ji,jj), pQsen(ji,jj), pQlat(ji,jj),             & 
     629            &                    pEvap=zevap, prhoa=zrho, pfact_evap=zfact_evap       ) 
     630 
     631         IF( PRESENT(pEvap) ) pEvap(ji,jj) = zevap 
    584632         IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho 
    585  
     633    
    586634      END_2D 
    587635   END SUBROUTINE BULK_FORMULA_VCTR 
    588  
    589  
    590    SUBROUTINE BULK_FORMULA_SCLR( pzu, pTs, pqs, pTa, pqa, & 
    591       &                          pCd, pCh, pCe,           & 
    592       &                          pwnd, pUb, pslp,         & 
    593       &                          pTau, pQsen, pQlat,  pEvap, prhoa ) 
    594       !!---------------------------------------------------------------------------------- 
    595       REAL(wp),                     INTENT(in)  :: pzu  ! height above the sea-level where all this takes place (normally 10m) 
    596       REAL(wp), INTENT(in)  :: pTs  ! water temperature at the air-sea interface [K] 
    597       REAL(wp), INTENT(in)  :: pqs  ! satur. spec. hum. at T=pTs   [kg/kg] 
    598       REAL(wp), INTENT(in)  :: pTa  ! potential air temperature at z=pzu [K] 
    599       REAL(wp), INTENT(in)  :: pqa  ! specific humidity at z=pzu [kg/kg] 
    600       REAL(wp), INTENT(in)  :: pCd 
    601       REAL(wp), INTENT(in)  :: pCh 
    602       REAL(wp), INTENT(in)  :: pCe 
    603       REAL(wp), INTENT(in)  :: pwnd ! wind speed module at z=pzu [m/s] 
    604       REAL(wp), INTENT(in)  :: pUb  ! bulk wind speed at z=pzu (inc. pot. effect of gustiness etc) [m/s] 
    605       REAL(wp), INTENT(in)  :: pslp ! sea-level atmospheric pressure [Pa] 
    606       !! 
    607       REAL(wp), INTENT(out) :: pTau  ! module of the wind stress [N/m^2] 
    608       REAL(wp), INTENT(out) :: pQsen !  [W/m^2] 
    609       REAL(wp), INTENT(out) :: pQlat !  [W/m^2] 
    610       !! 
    611       REAL(wp), INTENT(out), OPTIONAL :: pEvap ! Evaporation [kg/m^2/s] 
    612       REAL(wp), INTENT(out), OPTIONAL :: prhoa ! Air density at z=pzu [kg/m^3] 
    613       !! 
    614       REAL(wp) :: ztaa, zgamma, zrho, zUrho, zevap 
    615       INTEGER  :: jq 
    616       !!---------------------------------------------------------------------------------- 
    617  
    618       !! Need ztaa, absolute temperature at pzu (formula to estimate rho_air needs absolute temperature, not the potential temperature "pTa") 
    619       ztaa = pTa ! first guess... 
    620       DO jq = 1, 4 
    621          zgamma = gamma_moist( 0.5*(ztaa+pTs) , pqa ) 
    622          ztaa = pTa - zgamma*pzu   ! Absolute temp. is slightly colder... 
    623       END DO 
    624       zrho = rho_air(ztaa, pqa, pslp) 
    625       zrho = rho_air(ztaa, pqa, pslp-zrho*grav*pzu) ! taking into account that we are pzu m above the sea level where SLP is given! 
    626  
    627       zUrho = pUb*MAX(zrho, 1._wp)     ! rho*U10 
    628  
    629       pTau = zUrho * pCd * pwnd ! Wind stress module 
    630  
    631       zevap = zUrho * pCe * (pqa - pqs) 
    632       pQsen = zUrho * pCh * (pTa - pTs) * cp_air(pqa) 
    633       pQlat = L_vap(pTs) * zevap 
    634  
    635       IF( PRESENT(pEvap) ) pEvap = - zevap 
    636       IF( PRESENT(prhoa) ) prhoa = zrho 
    637  
    638    END SUBROUTINE BULK_FORMULA_SCLR 
    639  
    640  
    641636 
    642637 
    … …  
    645640      !!                           ***  FUNCTION alpha_sw_vctr  *** 
    646641      !! 
    647       !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa) 
     642      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (i.e. P =~ 101000 Pa) 
    648643      !! 
    649644      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 
    … …  
    659654      !!                           ***  FUNCTION alpha_sw_sclr  *** 
    660655      !! 
    661       !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (P =~ 1010 hpa) 
     656      !! ** Purpose : ROUGH estimate of the thermal expansion coefficient of sea-water at the surface (i.e. P =~ 101000 Pa) 
    662657      !! 
    663658      !! ** Author: L. Brodeau, june 2016 / AeroBulk (https://github.com/brodeau/aerobulk/) 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcblk_skin_coare.F90

    r12511 r13540  
    8989      REAL(wp) :: zQabs, zdlt, zfr, zalfa, zqlat, zus 
    9090      !!--------------------------------------------------------------------- 
    91       DO_2D_11_11 
     91      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    9292 
    9393         zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta, 
    … …  
    156156      ztime = REAL(nsec_day,wp)/(24._wp*3600._wp) ! time of current time step since 00:00 for current day (UTC) -> ztime = 0 -> 00:00 / ztime = 0.5 -> 12:00 ... 
    157157 
    158       DO_2D_11_11 
     158      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    159159 
    160160         l_exit       = .FALSE. 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcblk_skin_ecmwf.F90

    r12511 r13540  
    9595      REAL(wp) :: zQabs, zdlt, zfr, zalfa, zus 
    9696      !!--------------------------------------------------------------------- 
    97       DO_2D_11_11 
     97      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    9898 
    9999         zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta, 
    … …  
    173173      IF( PRESENT(pustk) ) l_pustk_known = .TRUE. 
    174174 
    175       DO_2D_11_11 
     175      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    176176 
    177177         zHwl = Hz_wl(ji,jj) ! first guess for warm-layer depth (and unique..., less advanced than COARE3p6 !) 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbccpl.F90

    r12511 r13540  
    4141#endif 
    4242#if defined key_si3 
    43    USE icethd_dh      ! for CALL ice_thd_snwblow 
     43   USE icevar         ! for CALL ice_var_snwblow 
    4444#endif 
    4545   ! 
    … …  
    4848   USE lib_mpp        ! distribued memory computing library 
    4949   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     50 
     51#if defined key_oasis3  
     52   USE mod_oasis, ONLY : OASIS_Sent, OASIS_ToRest, OASIS_SentOut, OASIS_ToRestOut  
     53#endif  
    5054 
    5155   IMPLICIT NONE 
    … …  
    152156   INTEGER, PARAMETER ::   jps_wlev   = 32   ! water level  
    153157   INTEGER, PARAMETER ::   jps_fice1  = 33   ! first-order ice concentration (for semi-implicit coupling of atmos-ice fluxes) 
    154    INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area 
     158   INTEGER, PARAMETER ::   jps_a_p    = 34   ! meltpond area fraction 
    155159   INTEGER, PARAMETER ::   jps_ht_p   = 35   ! meltpond thickness 
    156160   INTEGER, PARAMETER ::   jps_kice   = 36   ! sea ice effective conductivity 
    … …  
    159163 
    160164   INTEGER, PARAMETER ::   jpsnd      = 38   ! total number of fields sent  
     165 
     166#if ! defined key_oasis3  
     167   ! Dummy variables to enable compilation when oasis3 is not being used  
     168   INTEGER                    ::   OASIS_Sent        = -1  
     169   INTEGER                    ::   OASIS_SentOut     = -1  
     170   INTEGER                    ::   OASIS_ToRest      = -1  
     171   INTEGER                    ::   OASIS_ToRestOut   = -1  
     172#endif  
    161173 
    162174   !                                  !!** namelist namsbc_cpl ** 
    … …  
    184196   LOGICAL     ::   ln_usecplmask         !  use a coupling mask file to merge data received from several models 
    185197                                          !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 
     198   LOGICAL     ::   ln_scale_ice_flux     !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration)  
     199 
    186200   TYPE ::   DYNARR      
    187201      REAL(wp), POINTER, DIMENSION(:,:,:) ::   z3    
    … …  
    191205 
    192206   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky) 
     207#if defined key_si3 || defined key_cice 
     208   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i_last_couple !: Ice fractional area at last coupling time 
     209#endif 
    193210 
    194211   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]  
    … …  
    199216   !! Substitution 
    200217#  include "do_loop_substitute.h90" 
     218#  include "domzgr_substitute.h90" 
    201219   !!---------------------------------------------------------------------- 
    202220   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    210228      !!             ***  FUNCTION sbc_cpl_alloc  *** 
    211229      !!---------------------------------------------------------------------- 
    212       INTEGER :: ierr(4) 
     230      INTEGER :: ierr(5) 
    213231      !!---------------------------------------------------------------------- 
    214232      ierr(:) = 0 
    … …  
    220238#endif 
    221239      ALLOCATE( xcplmask(jpi,jpj,0:nn_cplmodel) , STAT=ierr(3) ) 
    222       ! 
    223       IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(4) )  
     240#if defined key_si3 || defined key_cice 
     241      ALLOCATE( a_i_last_couple(jpi,jpj,jpl) , STAT=ierr(4) ) 
     242#endif 
     243      ! 
     244      IF( .NOT. ln_apr_dyn ) ALLOCATE( ssh_ib(jpi,jpj), ssh_ibb(jpi,jpj), apr(jpi, jpj), STAT=ierr(5) ) 
    224245 
    225246      sbc_cpl_alloc = MAXVAL( ierr ) 
    … …  
    248269      REAL(wp), DIMENSION(jpi,jpj) ::   zacs, zaos 
    249270      !! 
    250       NAMELIST/namsbc_cpl/  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2  ,   &  
     271      NAMELIST/namsbc_cpl/  nn_cplmodel  , ln_usecplmask, nn_cats_cpl , ln_scale_ice_flux,             & 
     272         &                  sn_snd_temp  , sn_snd_alb   , sn_snd_thick, sn_snd_crt   , sn_snd_co2   ,  &  
    251273         &                  sn_snd_ttilyr, sn_snd_cond  , sn_snd_mpnd , sn_snd_sstfrz, sn_snd_thick1,  &  
    252          &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc,   &  
    253          &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr  ,   &  
     274         &                  sn_snd_ifrac , sn_snd_crtw  , sn_snd_wlev , sn_rcv_hsig  , sn_rcv_phioc ,  &  
     275         &                  sn_rcv_w10m  , sn_rcv_taumod, sn_rcv_tau  , sn_rcv_dqnsdt, sn_rcv_qsr   ,  &  
    254276         &                  sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum  , sn_rcv_tauwoc,  & 
    255          &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal  ,   & 
    256          &                  sn_rcv_iceflx, sn_rcv_co2   , nn_cplmodel , ln_usecplmask, sn_rcv_mslp ,   & 
    257          &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq , sn_rcv_tauw, nn_cats_cpl  ,   & 
     277         &                  sn_rcv_wdrag , sn_rcv_qns   , sn_rcv_emp  , sn_rcv_rnf   , sn_rcv_cal   ,  & 
     278         &                  sn_rcv_iceflx, sn_rcv_co2   , sn_rcv_mslp ,                                & 
     279         &                  sn_rcv_icb   , sn_rcv_isf   , sn_rcv_wfreq, sn_rcv_tauw  ,                 & 
    258280         &                  sn_rcv_ts_ice 
    259  
    260281      !!--------------------------------------------------------------------- 
    261282      ! 
    … …  
    277298      ENDIF 
    278299      IF( lwp .AND. ln_cpl ) THEN                        ! control print 
     300         WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
     301         WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
     302         WRITE(numout,*)'  ln_scale_ice_flux                   = ', ln_scale_ice_flux 
     303         WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl 
    279304         WRITE(numout,*)'  received fields (mutiple ice categogies)' 
    280305         WRITE(numout,*)'      10m wind module                 = ', TRIM(sn_rcv_w10m%cldes  ), ' (', TRIM(sn_rcv_w10m%clcat  ), ')' 
    … …  
    325350         WRITE(numout,*)'                      - orientation   = ', sn_snd_crtw%clvor  
    326351         WRITE(numout,*)'                      - mesh          = ', sn_snd_crtw%clvgrd  
    327          WRITE(numout,*)'  nn_cplmodel                         = ', nn_cplmodel 
    328          WRITE(numout,*)'  ln_usecplmask                       = ', ln_usecplmask 
    329          WRITE(numout,*)'  nn_cats_cpl                         = ', nn_cats_cpl 
    330352      ENDIF 
    331353 
    … …  
    364386      !  
    365387      ! Vectors: change of sign at north fold ONLY if on the local grid 
    366       IF( TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM(sn_rcv_tau%cldes ) == 'oce and ice') THEN ! avoid working with the atmospheric fields if they are not coupled 
     388      IF(       TRIM( sn_rcv_tau%cldes ) == 'oce only' .OR. TRIM( sn_rcv_tau%cldes ) == 'oce and ice'  & 
     389           .OR. TRIM( sn_rcv_tau%cldes ) == 'mixed oce-ice' ) THEN ! avoid working with the atmospheric fields if they are not coupled 
     390      ! 
    367391      IF( TRIM( sn_rcv_tau%clvor ) == 'local grid' )   srcv(jpr_otx1:jpr_itz2)%nsgn = -1. 
    368392       
    … …  
    695719         ! Change first letter to couple with atmosphere if already coupled OPA 
    696720         ! this is nedeed as each variable name used in the namcouple must be unique: 
    697          ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere 
     721         ! for example O_Runoff received by OPA from SAS and therefore S_Runoff received by SAS from the Atmosphere 
    698722         DO jn = 1, jprcv 
    699723            IF( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname)) 
    … …  
    819843      END SELECT 
    820844 
     845      ! Initialise ice fractions from last coupling time to zero (needed by Met-Office) 
     846#if defined key_si3 || defined key_cice 
     847       a_i_last_couple(:,:,:) = 0._wp 
     848#endif 
    821849      !                                                      ! ------------------------- !  
    822850      !                                                      !      Ice Meltponds        !  
    … …  
    10361064         xcplmask(:,:,:) = 0. 
    10371065         CALL iom_open( 'cplmask', inum ) 
    1038          CALL iom_get( inum, jpdom_unknown, 'cplmask', xcplmask(1:nlci,1:nlcj,1:nn_cplmodel),   & 
    1039             &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,nn_cplmodel /) ) 
     1066         CALL iom_get( inum, jpdom_unknown, 'cplmask', xcplmask(1:jpi,1:jpj,1:nn_cplmodel),   & 
     1067            &          kstart = (/ mig(1),mjg(1),1 /), kcount = (/ jpi,jpj,nn_cplmodel /) ) 
    10401068         CALL iom_close( inum ) 
    10411069      ELSE 
    … …  
    11071135      REAL(wp) ::   zcdrag = 1.5e-3        ! drag coefficient 
    11081136      REAL(wp) ::   zzx, zzy               ! temporary variables 
    1109       REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr 
     1137      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty, zmsk, zemp, zqns, zqsr, zcloud_fra 
    11101138      !!---------------------------------------------------------------------- 
    11111139      ! 
    … …  
    11151143         IF( ln_dm2dc .AND. ncpl_qsr_freq /= 86400 )   & 
    11161144            &   CALL ctl_stop( 'sbc_cpl_rcv: diurnal cycle reconstruction (ln_dm2dc) needs daily couping for solar radiation' ) 
    1117          ncpl_qsr_freq = 86400 / ncpl_qsr_freq   ! used by top 
     1145 
     1146         IF( ncpl_qsr_freq /= 0) ncpl_qsr_freq = 86400 / ncpl_qsr_freq ! used by top 
     1147          
    11181148      ENDIF 
    11191149      ! 
    … …  
    11651195            !                               
    11661196            IF( srcv(jpr_otx1)%clgrid == 'T' ) THEN 
    1167                DO_2D_00_00 
     1197               DO_2D( 0, 0, 0, 0 )                                        ! T ==> (U,V) 
    11681198                  frcv(jpr_otx1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_otx1)%z3(ji+1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) ) 
    11691199                  frcv(jpr_oty1)%z3(ji,jj,1) = 0.5 * ( frcv(jpr_oty1)%z3(ji  ,jj+1,1) + frcv(jpr_oty1)%z3(ji,jj,1) ) 
    11701200               END_2D 
    1171                CALL lbc_lnk_multi( 'sbccpl', frcv(jpr_otx1)%z3(:,:,1), 'U',  -1., frcv(jpr_oty1)%z3(:,:,1), 'V',  -1. ) 
     1201               CALL lbc_lnk_multi( 'sbccpl', frcv(jpr_otx1)%z3(:,:,1), 'U',  -1.0_wp, frcv(jpr_oty1)%z3(:,:,1), 'V',  -1.0_wp ) 
    11721202            ENDIF 
    11731203            llnewtx = .TRUE. 
    … …  
    11891219         ! => need to be done only when otx1 was changed 
    11901220         IF( llnewtx ) THEN 
    1191             DO_2D_00_00 
     1221            DO_2D( 0, 0, 0, 0 ) 
    11921222               zzx = frcv(jpr_otx1)%z3(ji-1,jj  ,1) + frcv(jpr_otx1)%z3(ji,jj,1) 
    11931223               zzy = frcv(jpr_oty1)%z3(ji  ,jj-1,1) + frcv(jpr_oty1)%z3(ji,jj,1) 
    11941224               frcv(jpr_taum)%z3(ji,jj,1) = 0.5 * SQRT( zzx * zzx + zzy * zzy ) 
    11951225            END_2D 
    1196             CALL lbc_lnk( 'sbccpl', frcv(jpr_taum)%z3(:,:,1), 'T', 1. ) 
     1226            CALL lbc_lnk( 'sbccpl', frcv(jpr_taum)%z3(:,:,1), 'T', 1.0_wp ) 
    11971227            llnewtau = .TRUE. 
    11981228         ELSE 
    … …  
    12141244         IF( llnewtau ) THEN  
    12151245            zcoef = 1. / ( zrhoa * zcdrag )  
    1216             DO_2D_11_11 
     1246            DO_2D( 1, 1, 1, 1 ) 
    12171247               frcv(jpr_w10m)%z3(ji,jj,1) = SQRT( frcv(jpr_taum)%z3(ji,jj,1) * zcoef ) 
    12181248            END_2D 
    12191249         ENDIF 
    12201250      ENDIF 
    1221  
     1251!!$      !                                                      ! ========================= ! 
     1252!!$      SELECT CASE( TRIM( sn_rcv_clouds%cldes ) )             !       cloud fraction      ! 
     1253!!$      !                                                      ! ========================= ! 
     1254!!$      cloud_fra(:,:) = frcv(jpr_clfra)*z3(:,:,1) 
     1255!!$      END SELECT 
     1256!!$ 
     1257      zcloud_fra(:,:) = pp_cldf   ! should be real cloud fraction instead (as in the bulk) but needs to be read from atm. 
     1258      IF( ln_mixcpl ) THEN 
     1259         cloud_fra(:,:) = cloud_fra(:,:) * xcplmask(:,:,0) + zcloud_fra(:,:)* zmsk(:,:) 
     1260      ELSE 
     1261         cloud_fra(:,:) = zcloud_fra(:,:) 
     1262      ENDIF 
     1263      !                                                      ! ========================= ! 
    12221264      ! u(v)tau and taum will be modified by ice model 
    12231265      ! -> need to be reset before each call of the ice/fsbc       
    … …  
    14791521      INTEGER ::   ji, jj   ! dummy loop indices 
    14801522      INTEGER ::   itx      ! index of taux over ice 
     1523      REAL(wp)                     ::   zztmp1, zztmp2 
    14811524      REAL(wp), DIMENSION(jpi,jpj) ::   ztx, zty  
    14821525      !!---------------------------------------------------------------------- 
    … …  
    15421585            p_taui(:,:) = frcv(jpr_itx1)%z3(:,:,1)                   ! (U,V) ==> (U,V) 
    15431586            p_tauj(:,:) = frcv(jpr_ity1)%z3(:,:,1) 
    1544          CASE( 'F' ) 
    1545             DO_2D_00_00 
    1546                p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji,jj,1) + frcv(jpr_itx1)%z3(ji  ,jj-1,1) ) 
    1547                p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji,jj,1) + frcv(jpr_ity1)%z3(ji-1,jj  ,1) ) 
     1587         CASE( 'T' ) 
     1588            DO_2D( 0, 0, 0, 0 )                    ! T ==> (U,V) 
     1589               ! take care of the land-sea mask to avoid "pollution" of coastal stress. p[uv]taui used in frazil and  rheology  
     1590               zztmp1 = 0.5_wp * ( 2. - umask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji+1,jj  ,1) ) 
     1591               zztmp2 = 0.5_wp * ( 2. - vmask(ji,jj,1) ) * MAX( tmask(ji,jj,1),tmask(ji  ,jj+1,1) ) 
     1592               p_taui(ji,jj) = zztmp1 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) ) 
     1593               p_tauj(ji,jj) = zztmp2 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) ) 
    15481594            END_2D 
    1549          CASE( 'T' ) 
    1550             DO_2D_00_00 
    1551                p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj  ,1) + frcv(jpr_itx1)%z3(ji,jj,1) ) 
    1552                p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji  ,jj+1,1) + frcv(jpr_ity1)%z3(ji,jj,1) ) 
    1553             END_2D 
    1554          CASE( 'I' ) 
    1555             DO_2D_00_00 
    1556                p_taui(ji,jj) = 0.5 * ( frcv(jpr_itx1)%z3(ji+1,jj+1,1) + frcv(jpr_itx1)%z3(ji+1,jj  ,1) ) 
    1557                p_tauj(ji,jj) = 0.5 * ( frcv(jpr_ity1)%z3(ji+1,jj+1,1) + frcv(jpr_ity1)%z3(ji  ,jj+1,1) ) 
    1558             END_2D 
     1595            CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. ) 
    15591596         END SELECT 
    1560          IF( srcv(jpr_itx1)%clgrid /= 'U' ) THEN  
    1561             CALL lbc_lnk_multi( 'sbccpl', p_taui, 'U',  -1., p_tauj, 'V',  -1. ) 
    1562          ENDIF 
    15631597          
    15641598      ENDIF 
    … …  
    16261660      ! 
    16271661      INTEGER  ::   ji, jj, jl   ! dummy loop index 
    1628       REAL(wp) ::   ztri         ! local scalar 
    16291662      REAL(wp), DIMENSION(jpi,jpj)     ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw 
    16301663      REAL(wp), DIMENSION(jpi,jpj)     ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip  , zevap_oce, zdevap_ice 
    16311664      REAL(wp), DIMENSION(jpi,jpj)     ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
     1665      REAL(wp), DIMENSION(jpi,jpj)     ::   zevap_ice_total 
    16321666      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zevap_ice, zqtr_ice_top, ztsu 
     1667      REAL(wp), DIMENSION(jpi,jpj)     ::   ztri 
    16331668      !!---------------------------------------------------------------------- 
    16341669      ! 
    … …  
    16501685         ztprecip(:,:) =   frcv(jpr_rain)%z3(:,:,1) + zsprecip(:,:)  ! May need to ensure positive here 
    16511686         zemp_tot(:,:) =   frcv(jpr_tevp)%z3(:,:,1) - ztprecip(:,:) 
    1652          zemp_ice(:,:) = ( frcv(jpr_ievp)%z3(:,:,1) - frcv(jpr_snow)%z3(:,:,1) ) * picefr(:,:) 
    16531687      CASE( 'oce and ice'   )   ! received fields: jpr_sbpr, jpr_semp, jpr_oemp, jpr_ievp 
    16541688         zemp_tot(:,:) = ziceld(:,:) * frcv(jpr_oemp)%z3(:,:,1) + picefr(:,:) * frcv(jpr_sbpr)%z3(:,:,1) 
    … …  
    16621696 
    16631697#if defined key_si3 
     1698 
     1699      ! --- evaporation over ice (kg/m2/s) --- ! 
     1700      IF (ln_scale_ice_flux) THEN ! typically met-office requirements 
     1701         IF (sn_rcv_emp%clcat == 'yes') THEN 
     1702            WHERE( a_i(:,:,:) > 1.e-10 )  ; zevap_ice(:,:,:) = frcv(jpr_ievp)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     1703            ELSEWHERE                     ; zevap_ice(:,:,:) = 0._wp 
     1704            END WHERE 
     1705            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:) 
     1706            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp 
     1707            END WHERE 
     1708         ELSE 
     1709            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) * SUM( a_i_last_couple, dim=3 ) / picefr(:,:) 
     1710            ELSEWHERE                     ; zevap_ice(:,:,1) = 0._wp 
     1711            END WHERE 
     1712            zevap_ice_total(:,:) = zevap_ice(:,:,1) 
     1713            DO jl = 2, jpl 
     1714               zevap_ice(:,:,jl) = zevap_ice(:,:,1) 
     1715            ENDDO 
     1716         ENDIF 
     1717      ELSE 
     1718         IF (sn_rcv_emp%clcat == 'yes') THEN 
     1719            zevap_ice(:,:,1:jpl) = frcv(jpr_ievp)%z3(:,:,1:jpl) 
     1720            WHERE( picefr(:,:) > 1.e-10 ) ; zevap_ice_total(:,:) = SUM( zevap_ice(:,:,:) * a_i(:,:,:), dim=3 ) / picefr(:,:) 
     1721            ELSEWHERE                     ; zevap_ice_total(:,:) = 0._wp 
     1722            END WHERE 
     1723         ELSE 
     1724            zevap_ice(:,:,1) = frcv(jpr_ievp)%z3(:,:,1) 
     1725            zevap_ice_total(:,:) = zevap_ice(:,:,1) 
     1726            DO jl = 2, jpl 
     1727               zevap_ice(:,:,jl) = zevap_ice(:,:,1) 
     1728            ENDDO 
     1729         ENDIF 
     1730      ENDIF 
     1731 
     1732      IF ( TRIM( sn_rcv_emp%cldes ) == 'conservative' ) THEN 
     1733         ! For conservative case zemp_ice has not been defined yet. Do it now. 
     1734         zemp_ice(:,:) = zevap_ice_total(:,:) * picefr(:,:) - frcv(jpr_snow)%z3(:,:,1) * picefr(:,:) 
     1735      ENDIF 
     1736 
    16641737      ! zsnw = snow fraction over ice after wind blowing (=picefr if no blowing) 
    1665       zsnw(:,:) = 0._wp   ;   CALL ice_thd_snwblow( ziceld, zsnw ) 
     1738      zsnw(:,:) = 0._wp   ;   CALL ice_var_snwblow( ziceld, zsnw ) 
    16661739       
    16671740      ! --- evaporation minus precipitation corrected (because of wind blowing on snow) --- ! 
    … …  
    16701743 
    16711744      ! --- evaporation over ocean (used later for qemp) --- ! 
    1672       zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) 
    1673  
    1674       ! --- evaporation over ice (kg/m2/s) --- ! 
    1675       DO jl=1,jpl 
    1676          IF(sn_rcv_emp%clcat == 'yes') THEN   ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl) 
    1677          ELSE                                  ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,1 )   ;   ENDIF 
    1678       ENDDO 
     1745      zevap_oce(:,:) = frcv(jpr_tevp)%z3(:,:,1) - zevap_ice_total(:,:) * picefr(:,:) 
    16791746 
    16801747      ! since the sensitivity of evap to temperature (devap/dT) is not prescribed by the atmosphere, we set it to 0 
    … …  
    17541821!!      IF( srcv(jpr_rnf)%laction )   CALL iom_put( 'runoffs' , rnf(:,:) * tmask(:,:,1)                                 )  ! runoff 
    17551822!!      IF( srcv(jpr_isf)%laction )   CALL iom_put( 'iceshelf_cea', -fwfisf(:,:) * tmask(:,:,1)                         )  ! iceshelf 
    1756       IF( srcv(jpr_cal)%laction )   CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving 
    1757       IF( srcv(jpr_icb)%laction )   CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs 
    1758       IF( iom_use('snowpre') )      CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow 
    1759       IF( iom_use('precip') )       CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation 
    1760       IF( iom_use('rain') )         CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation  
    1761       IF( iom_use('snow_ao_cea') )  CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average) 
    1762       IF( iom_use('snow_ai_cea') )  CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average) 
    1763       IF( iom_use('rain_ao_cea') )  CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average) 
    1764       IF( iom_use('subl_ai_cea') )  CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )  ! Sublimation over sea-ice (cell average) 
    1765       IF( iom_use('evap_ao_cea') )  CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  & 
    1766          &                                                        - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average) 
     1823      IF( srcv(jpr_cal)%laction )    CALL iom_put( 'calving_cea' , frcv(jpr_cal)%z3(:,:,1) * tmask(:,:,1)                )  ! calving 
     1824      IF( srcv(jpr_icb)%laction )    CALL iom_put( 'iceberg_cea' , frcv(jpr_icb)%z3(:,:,1) * tmask(:,:,1)                )  ! icebergs 
     1825      IF( iom_use('snowpre') )       CALL iom_put( 'snowpre'     , sprecip(:,:)                                          )  ! Snow 
     1826      IF( iom_use('precip') )        CALL iom_put( 'precip'      , tprecip(:,:)                                          )  ! total  precipitation 
     1827      IF( iom_use('rain') )          CALL iom_put( 'rain'        , tprecip(:,:) - sprecip(:,:)                           )  ! liquid precipitation  
     1828      IF( iom_use('snow_ao_cea') )   CALL iom_put( 'snow_ao_cea' , sprecip(:,:) * ( 1._wp - zsnw(:,:) )                  )  ! Snow over ice-free ocean  (cell average) 
     1829      IF( iom_use('snow_ai_cea') )   CALL iom_put( 'snow_ai_cea' , sprecip(:,:) *           zsnw(:,:)                    )  ! Snow over sea-ice         (cell average) 
     1830      IF( iom_use('rain_ao_cea') )   CALL iom_put( 'rain_ao_cea' , ( tprecip(:,:) - sprecip(:,:) ) * picefr(:,:)         )  ! liquid precipitation over ocean (cell average) 
     1831      IF( iom_use('subl_ai_cea') )   CALL iom_put( 'subl_ai_cea' , frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) * tmask(:,:,1) )     ! Sublimation over sea-ice (cell average) 
     1832      IF( iom_use('evap_ao_cea') )   CALL iom_put( 'evap_ao_cea' , ( frcv(jpr_tevp)%z3(:,:,1)  & 
     1833         &                                                         - frcv(jpr_ievp)%z3(:,:,1) * picefr(:,:) ) * tmask(:,:,1) ) ! ice-free oce evap (cell average) 
    17671834      ! note: runoff output is done in sbcrnf (which includes icebergs too) and iceshelf output is done in sbcisf 
    17681835      ! 
    … …  
    17721839      CASE( 'oce only' )         ! the required field is directly provided 
    17731840         zqns_tot(:,:) = frcv(jpr_qnsoce)%z3(:,:,1) 
     1841         ! For Met Office sea ice non-solar fluxes are already delt with by JULES so setting to zero 
     1842         ! here so the only flux is the ocean only one. 
     1843         zqns_ice(:,:,:) = 0._wp  
    17741844      CASE( 'conservative' )     ! the required fields are directly provided 
    17751845         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1) 
    … …  
    17891859            ENDDO 
    17901860         ELSE 
    1791             qns_tot(:,:) = qns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1) 
     1861            zqns_tot(:,:) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1) 
    17921862            DO jl = 1, jpl 
    1793                zqns_tot(:,:   ) = zqns_tot(:,:) + picefr(:,:) * frcv(jpr_qnsice)%z3(:,:,1) 
    17941863               zqns_ice(:,:,jl) = frcv(jpr_qnsice)%z3(:,:,1) 
    17951864            END DO 
    … …  
    18021871               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:,jl)    & 
    18031872                  &             + frcv(jpr_dqnsdt)%z3(:,:,jl) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   & 
    1804                   &                                                                + pist(:,:,jl) * picefr(:,:) ) ) 
     1873                  &                                             + pist(:,:,jl) * picefr(:,:) ) ) 
    18051874            END DO 
    18061875         ELSE 
    … …  
    18081877               zqns_ice(:,:,jl) = frcv(jpr_qnsmix)%z3(:,:, 1)    & 
    18091878                  &             + frcv(jpr_dqnsdt)%z3(:,:, 1) * ( pist(:,:,jl) - ( ( rt0 + psst(:,:) ) * ziceld(:,:)   & 
    1810                   &                                                                + pist(:,:,jl) * picefr(:,:) ) ) 
     1879                  &                                             + pist(:,:,jl) * picefr(:,:) ) ) 
    18111880            END DO 
    18121881         ENDIF 
    … …  
    19141983      CASE( 'oce only' ) 
    19151984         zqsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) ) 
     1985         ! For Met Office sea ice solar fluxes are already delt with by JULES so setting to zero 
     1986         ! here so the only flux is the ocean only one. 
     1987         zqsr_ice(:,:,:) = 0._wp 
    19161988      CASE( 'conservative' ) 
    19171989         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1) 
    … …  
    19322004            END DO 
    19332005         ELSE 
    1934             qsr_tot(:,:   ) = qsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1) 
     2006            zqsr_tot(:,:) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1) 
    19352007            DO jl = 1, jpl 
    1936                zqsr_tot(:,:   ) = zqsr_tot(:,:) + picefr(:,:) * frcv(jpr_qsrice)%z3(:,:,1) 
    19372008               zqsr_ice(:,:,jl) = frcv(jpr_qsrice)%z3(:,:,1) 
    19382009            END DO 
    … …  
    20002071            ENDDO 
    20012072         ENDIF 
     2073      CASE( 'none' )  
     2074         zdqns_ice(:,:,:) = 0._wp 
    20022075      END SELECT 
    20032076       
    … …  
    20152088      !                                                      ! ========================= ! 
    20162089      CASE ('coupled') 
    2017          IF( ln_mixcpl ) THEN 
    2018             DO jl=1,jpl 
    2019                qml_ice(:,:,jl) = qml_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_topm)%z3(:,:,jl) * zmsk(:,:) 
    2020                qcn_ice(:,:,jl) = qcn_ice(:,:,jl) * xcplmask(:,:,0) + frcv(jpr_botm)%z3(:,:,jl) * zmsk(:,:) 
    2021             ENDDO 
     2090         IF (ln_scale_ice_flux) THEN 
     2091            WHERE( a_i(:,:,:) > 1.e-10_wp ) 
     2092               qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     2093               qcn_ice(:,:,:) = frcv(jpr_botm)%z3(:,:,:) * a_i_last_couple(:,:,:) / a_i(:,:,:) 
     2094            ELSEWHERE 
     2095               qml_ice(:,:,:) = 0.0_wp 
     2096               qcn_ice(:,:,:) = 0.0_wp 
     2097            END WHERE 
    20222098         ELSE 
    20232099            qml_ice(:,:,:) = frcv(jpr_topm)%z3(:,:,:) 
    … …  
    20302106      IF( .NOT.ln_cndflx ) THEN                              !==  No conduction flux as surface forcing  ==! 
    20312107         ! 
    2032          !                    ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
    2033          ztri = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice    ! surface transmission when hi>10cm (Grenfell Maykut 77) 
    2034          ! 
    2035          WHERE    ( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
    2036             zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ( ztri + ( 1._wp - ztri ) * ( 1._wp - phi(:,:,:) * 10._wp ) ) 
    2037          ELSEWHERE( phs(:,:,:) <= 0._wp .AND. phi(:,:,:) >= 0.1_wp )       ! constant (ztri) when hi>10cm 
    2038             zqtr_ice_top(:,:,:) = qsr_ice(:,:,:) * ztri 
    2039          ELSEWHERE                                                         ! zero when hs>0 
    2040             zqtr_ice_top(:,:,:) = 0._wp 
    2041          END WHERE 
     2108         IF( nn_qtrice == 0 ) THEN 
     2109            ! formulation derived from Grenfell and Maykut (1977), where transmission rate 
     2110            !    1) depends on cloudiness 
     2111            !       ! ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
     2112            !       !      should be real cloud fraction instead (as in the bulk) but needs to be read from atm. 
     2113            !    2) is 0 when there is any snow 
     2114            !    3) tends to 1 for thin ice 
     2115            ztri(:,:) = 0.18 * ( 1.0 - cloud_fra(:,:) ) + 0.35 * cloud_fra(:,:)  ! surface transmission when hi>10cm 
     2116            DO jl = 1, jpl 
     2117               WHERE    ( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) <  0.1_wp )       ! linear decrease from hi=0 to 10cm   
     2118                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ( ztri(:,:) + ( 1._wp - ztri(:,:) ) * ( 1._wp - phi(:,:,jl) * 10._wp ) ) 
     2119               ELSEWHERE( phs(:,:,jl) <= 0._wp .AND. phi(:,:,jl) >= 0.1_wp )       ! constant (ztri) when hi>10cm 
     2120                  zqtr_ice_top(:,:,jl) = zqsr_ice(:,:,jl) * ztri(:,:) 
     2121               ELSEWHERE                                                           ! zero when hs>0 
     2122                  zqtr_ice_top(:,:,jl) = 0._wp  
     2123               END WHERE 
     2124            ENDDO 
     2125         ELSEIF( nn_qtrice == 1 ) THEN 
     2126            ! formulation is derived from the thesis of M. Lebrun (2019). 
     2127            !    It represents the best fit using several sets of observations 
     2128            !    It comes with snow conductivities adapted to freezing/melting conditions (see icethd_zdf_bl99.F90) 
     2129            zqtr_ice_top(:,:,:) = 0.3_wp * zqsr_ice(:,:,:) 
     2130         ENDIF 
    20422131         !      
    20432132      ELSEIF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN      !==  conduction flux as surface forcing  ==! 
    20442133         ! 
    2045          !                    ! ===> here we must receive the qtr_ice_top array from the coupler 
    2046          !                           for now just assume zero (fully opaque ice) 
     2134         !          ! ===> here we must receive the qtr_ice_top array from the coupler 
     2135         !                 for now just assume zero (fully opaque ice) 
    20472136         zqtr_ice_top(:,:,:) = 0._wp 
    20482137         ! 
    … …  
    21012190      ! 
    21022191      isec = ( kt - nit000 ) * NINT( rn_Dt )        ! date of exchanges 
     2192      info = OASIS_idle 
    21032193 
    21042194      zfr_l(:,:) = 1.- fr_i(:,:) 
    … …  
    22392329      ENDIF 
    22402330 
     2331#if defined key_si3 || defined key_cice 
     2332      ! If this coupling was successful then save ice fraction for use between coupling points.  
     2333      ! This is needed for some calculations where the ice fraction at the last coupling point  
     2334      ! is needed.  
     2335      IF(  info == OASIS_Sent    .OR. info == OASIS_ToRest .OR. &  
     2336         & info == OASIS_SentOut .OR. info == OASIS_ToRestOut ) THEN  
     2337         IF ( sn_snd_thick%clcat == 'yes' ) THEN  
     2338           a_i_last_couple(:,:,1:jpl) = a_i(:,:,1:jpl) 
     2339         ENDIF 
     2340      ENDIF 
     2341#endif 
     2342 
    22412343      IF( ssnd(jps_fice1)%laction ) THEN 
    22422344         SELECT CASE( sn_snd_thick1%clcat ) 
    … …  
    23022404            SELECT CASE( sn_snd_mpnd%clcat )   
    23032405            CASE( 'yes' )   
    2304                ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl) 
     2406               ztmp3(:,:,1:jpl) =  a_ip_eff(:,:,1:jpl) 
    23052407               ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)   
    23062408            CASE( 'no' )   
    … …  
    23082410               ztmp4(:,:,:) = 0.0   
    23092411               DO jl=1,jpl   
    2310                  ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl)   
    2311                  ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl)  
     2412                 ztmp3(:,:,1) = ztmp3(:,:,1) + a_ip_frac(:,:,jpl) 
     2413                 ztmp4(:,:,1) = ztmp4(:,:,1) + h_ip(:,:,jpl) 
    23122414               ENDDO   
    23132415            CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_mpnd%clcat' )   
    … …  
    23702472            SELECT CASE( TRIM( sn_snd_crt%cldes ) ) 
    23712473            CASE( 'oce only'             )      ! C-grid ==> T 
    2372                DO_2D_00_00 
     2474               DO_2D( 0, 0, 0, 0 ) 
    23732475                  zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj  ,1,Kmm) ) 
    23742476                  zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji  ,jj-1,1,Kmm) )  
    23752477               END_2D 
    23762478            CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T   
    2377                DO_2D_00_00 
     2479               DO_2D( 0, 0, 0, 0 ) 
    23782480                  zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   
    23792481                  zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj) 
    … …  
    23812483                  zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  )     + v_ice(ji  ,jj-1  )     ) *  fr_i(ji,jj) 
    23822484               END_2D 
    2383                CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1., zity1, 'T', -1. ) 
     2485               CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1.0_wp, zity1, 'T', -1.0_wp ) 
    23842486            CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T 
    2385                DO_2D_00_00 
     2487               DO_2D( 0, 0, 0, 0 ) 
    23862488                  zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   & 
    23872489                     &         + 0.5 * ( u_ice(ji,jj  )     + u_ice(ji-1,jj    )     ) *  fr_i(ji,jj) 
    … …  
    23902492               END_2D 
    23912493            END SELECT 
    2392             CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.,  zoty1, ssnd(jps_ocy1)%clgrid, -1. ) 
     2494            CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocx1)%clgrid, -1.0_wp,  zoty1, ssnd(jps_ocy1)%clgrid, -1.0_wp ) 
    23932495            ! 
    23942496         ENDIF 
    … …  
    24472549          SELECT CASE( TRIM( sn_snd_crtw%cldes ) )  
    24482550          CASE( 'oce only'             )      ! C-grid ==> T  
    2449              DO_2D_00_00 
     2551             DO_2D( 0, 0, 0, 0 ) 
    24502552                zotx1(ji,jj) = 0.5 * ( uu(ji,jj,1,Kmm) + uu(ji-1,jj  ,1,Kmm) )  
    24512553                zoty1(ji,jj) = 0.5 * ( vv(ji,jj,1,Kmm) + vv(ji , jj-1,1,Kmm) )   
    24522554             END_2D 
    24532555          CASE( 'weighted oce and ice' )      ! Ocean and Ice on C-grid ==> T    
    2454              DO_2D_00_00 
     2556             DO_2D( 0, 0, 0, 0 ) 
    24552557                zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)    
    24562558                zoty1(ji,jj) = 0.5 * ( vv   (ji,jj,1,Kmm) + vv   (ji  ,jj-1,1,Kmm) ) * zfr_l(ji,jj)  
    … …  
    24582560                zity1(ji,jj) = 0.5 * ( v_ice(ji,jj  ) + v_ice(ji  ,jj-1  ) ) *  fr_i(ji,jj)  
    24592561             END_2D 
    2460              CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1.,  zity1, 'T', -1. )  
     2562             CALL lbc_lnk_multi( 'sbccpl', zitx1, 'T', -1.0_wp,  zity1, 'T', -1.0_wp )  
    24612563          CASE( 'mixed oce-ice'        )      ! Ocean and Ice on C-grid ==> T   
    2462              DO_2D_00_00 
     2564             DO_2D( 0, 0, 0, 0 ) 
    24632565                zotx1(ji,jj) = 0.5 * ( uu   (ji,jj,1,Kmm) + uu   (ji-1,jj  ,1,Kmm) ) * zfr_l(ji,jj)   &  
    24642566                   &         + 0.5 * ( u_ice(ji,jj  ) + u_ice(ji-1,jj    ) ) *  fr_i(ji,jj)  
    … …  
    24672569             END_2D 
    24682570          END SELECT 
    2469          CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocxw)%clgrid, -1., zoty1, ssnd(jps_ocyw)%clgrid, -1. )  
     2571         CALL lbc_lnk_multi( 'sbccpl', zotx1, ssnd(jps_ocxw)%clgrid, -1.0_wp, zoty1, ssnd(jps_ocyw)%clgrid, -1.0_wp )  
    24702572         !  
    24712573         !  

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcdcy.F90

    r12511 r13540  
    110110 
    111111      imask_night(:,:) = 0 
    112       DO_2D_11_11 
     112      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    113113         ztmpm = 0._wp 
    114114         IF( ABS(rab(ji,jj)) < 1. ) THEN         ! day duration is less than 24h 
    … …  
    193193 
    194194         zsin = SIN( zdecrad )   ;   zcos = COS( zdecrad ) 
    195          DO_2D_11_11 
     195         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    196196            ztmp = rad * gphit(ji,jj) 
    197197            raa(ji,jj) = SIN( ztmp ) * zsin 
    … …  
    202202         ! rab to test if the day time is equal to 0, less than 24h of full day 
    203203         rab(:,:) = -raa(:,:) / rbb(:,:) 
    204          DO_2D_11_11 
     204         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    205205            IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
    206206               ! When is it night? 
    … …  
    226226         !         Avoid possible infinite scaling factor, associated with very short daylight 
    227227         !         periods, by ignoring periods less than 1/1000th of a day (ticket #1040) 
    228          DO_2D_11_11 
     228         DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
    229229            IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
    230230               rscal(ji,jj) = 0.0_wp 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcflx.F90

    r12377 r13540  
    2929   PUBLIC sbc_flx       ! routine called by step.F90 
    3030 
    31    INTEGER , PARAMETER ::   jpfld   = 5   ! maximum number of files to read  
    3231   INTEGER , PARAMETER ::   jp_utau = 1   ! index of wind stress (i-component) file 
    3332   INTEGER , PARAMETER ::   jp_vtau = 2   ! index of wind stress (j-component) file 
    … …  
    3534   INTEGER , PARAMETER ::   jp_qsr  = 4   ! index of solar heat file 
    3635   INTEGER , PARAMETER ::   jp_emp  = 5   ! index of evaporation-precipation file 
     36 !!INTEGER , PARAMETER ::   jp_sfx  = 6   ! index of salt flux flux 
     37   INTEGER , PARAMETER ::   jpfld   = 5 !! 6 ! maximum number of files to read  
    3738   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf    ! structure of input fields (file informations, fields read) 
    3839 
    … …  
    5960      !!                   net downward radiative flux            qsr   (watt/m2) 
    6061      !!                   net upward freshwater (evapo - precip) emp   (kg/m2/s) 
     62      !!                   salt flux                              sfx   (pss*dh*rho/dt => g/m2/s) 
    6163      !! 
    6264      !!      CAUTION :  - never mask the surface stress fields 
    … …  
    7173      !!              - emp         upward mass flux (evap. - precip.) 
    7274      !!              - sfx         salt flux; set to zero at nit000 but possibly non-zero 
    73       !!                            if ice is present 
     75      !!                            if ice 
    7476      !!---------------------------------------------------------------------- 
    7577      INTEGER, INTENT(in) ::   kt   ! ocean time step 
    … …  
    8587      CHARACTER(len=100) ::  cn_dir                               ! Root directory for location of flx files 
    8688      TYPE(FLD_N), DIMENSION(jpfld) ::   slf_i                    ! array of namelist information structures 
    87       TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp ! informations about the fields to be read 
    88       NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp 
     89      TYPE(FLD_N) ::   sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp !!, sn_sfx ! informations about the fields to be read 
     90      NAMELIST/namsbc_flx/ cn_dir, sn_utau, sn_vtau, sn_qtot, sn_qsr, sn_emp !!, sn_sfx 
    8991      !!--------------------------------------------------------------------- 
    9092      ! 
    … …  
    105107         slf_i(jp_utau) = sn_utau   ;   slf_i(jp_vtau) = sn_vtau 
    106108         slf_i(jp_qtot) = sn_qtot   ;   slf_i(jp_qsr ) = sn_qsr  
    107          slf_i(jp_emp ) = sn_emp 
     109         slf_i(jp_emp ) = sn_emp !! ;   slf_i(jp_sfx ) = sn_sfx 
    108110         ! 
    109111         ALLOCATE( sf(jpfld), STAT=ierror )        ! set sf structure 
    … …  
    118120         CALL fld_fill( sf, slf_i, cn_dir, 'sbc_flx', 'flux formulation for ocean surface boundary condition', 'namsbc_flx' ) 
    119121         ! 
    120          sfx(:,:) = 0.0_wp                         ! salt flux due to freezing/melting (non-zero only if ice is present) 
    121          ! 
    122122      ENDIF 
    123123 
    … …  
    126126      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN                        ! update ocean fluxes at each SBC frequency 
    127127 
    128          IF( ln_dm2dc ) THEN   ;   qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) )   ! modify now Qsr to include the diurnal cycle 
    129          ELSE                  ;   qsr(:,:) =          sf(jp_qsr)%fnow(:,:,1) 
     128         IF( ln_dm2dc ) THEN   ! modify now Qsr to include the diurnal cycle 
     129            qsr(:,:) = sbc_dcy( sf(jp_qsr)%fnow(:,:,1) ) * tmask(ji,jj,1) 
     130         ELSE 
     131            DO_2D( 0, 0, 0, 0 ) 
     132               qsr(ji,jj) =          sf(jp_qsr)%fnow(ji,jj,1)   * tmask(ji,jj,1) 
     133            END_2D 
    130134         ENDIF 
    131          DO_2D_11_11 
    132             utau(ji,jj) = sf(jp_utau)%fnow(ji,jj,1) 
    133             vtau(ji,jj) = sf(jp_vtau)%fnow(ji,jj,1) 
    134             qns (ji,jj) = sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) 
    135             emp (ji,jj) = sf(jp_emp )%fnow(ji,jj,1) 
     135         DO_2D( 0, 0, 0, 0 )                                      ! set the ocean fluxes from read fields 
     136            utau(ji,jj) =   sf(jp_utau)%fnow(ji,jj,1)                              * umask(ji,jj,1) 
     137            vtau(ji,jj) =   sf(jp_vtau)%fnow(ji,jj,1)                              * vmask(ji,jj,1) 
     138            qns (ji,jj) = ( sf(jp_qtot)%fnow(ji,jj,1) - sf(jp_qsr)%fnow(ji,jj,1) ) * tmask(ji,jj,1) 
     139            emp (ji,jj) =   sf(jp_emp )%fnow(ji,jj,1)                              * tmask(ji,jj,1) 
     140            !!sfx (ji,jj) = sf(jp_sfx )%fnow(ji,jj,1)                              * tmask(ji,jj,1)  
    136141         END_2D 
    137142         !                                                        ! add to qns the heat due to e-p 
    138          qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST 
     143         !!clem: I do not think it is needed 
     144         !!qns(:,:) = qns(:,:) - emp(:,:) * sst_m(:,:) * rcp        ! mass flux is at SST 
    139145         ! 
    140          qns(:,:) = qns(:,:) * tmask(:,:,1) 
    141          emp(:,:) = emp(:,:) * tmask(:,:,1) 
     146         ! clem: without these lbc calls, it seems that the northfold is not ok (true in 3.6, not sure in 4.x)  
     147         CALL lbc_lnk_multi( 'sbcflx', utau, 'U', -1._wp, vtau, 'V', -1._wp, & 
     148            &                           qns, 'T',  1._wp, emp , 'T',  1._wp, qsr, 'T', 1._wp ) !! sfx, 'T', 1._wp  ) 
    142149         ! 
    143          !                                                        ! module of wind stress and wind speed at T-point 
    144          zcoef = 1. / ( zrhoa * zcdrag ) 
    145          DO_2D_00_00 
    146             ztx = utau(ji-1,jj  ) + utau(ji,jj)  
    147             zty = vtau(ji  ,jj-1) + vtau(ji,jj)  
    148             zmod = 0.5 * SQRT( ztx * ztx + zty * zty ) 
    149             taum(ji,jj) = zmod 
    150             wndm(ji,jj) = SQRT( zmod * zcoef ) 
    151          END_2D 
    152          taum(:,:) = taum(:,:) * tmask(:,:,1) ; wndm(:,:) = wndm(:,:) * tmask(:,:,1) 
    153          CALL lbc_lnk( 'sbcflx', taum(:,:), 'T', 1. )   ;   CALL lbc_lnk( 'sbcflx', wndm(:,:), 'T', 1. ) 
    154  
    155150         IF( nitend-nit000 <= 100 .AND. lwp ) THEN                ! control print (if less than 100 time-step asked) 
    156151            WRITE(numout,*)  
    … …  
    166161         ! 
    167162      ENDIF 
     163      !                                                           ! module of wind stress and wind speed at T-point 
     164      ! Note the use of 0.5*(2-umask) in order to unmask the stress along coastlines 
     165      zcoef = 1. / ( zrhoa * zcdrag ) 
     166      DO_2D( 0, 0, 0, 0 ) 
     167         ztx = ( utau(ji-1,jj  ) + utau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( umask(ji-1,jj  ,1), umask(ji,jj,1) ) ) 
     168         zty = ( vtau(ji  ,jj-1) + vtau(ji,jj) ) * 0.5_wp * ( 2._wp - MIN( vmask(ji  ,jj-1,1), vmask(ji,jj,1) ) )  
     169         zmod = 0.5_wp * SQRT( ztx * ztx + zty * zty ) * tmask(ji,jj,1) 
     170         taum(ji,jj) = zmod 
     171         wndm(ji,jj) = SQRT( zmod * zcoef )  !!clem: not used? 
     172      END_2D 
     173      ! 
     174      CALL lbc_lnk_multi( 'sbcflx', taum, 'T', 1._wp, wndm, 'T', 1._wp ) 
    168175      ! 
    169176   END SUBROUTINE sbc_flx 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcfwb.F90

    r12511 r13540  
    7171      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   ztmsk_tospread, zerp_cor    !   -      - 
    7272      REAL(wp)   ,DIMENSION(1) ::   z_fwfprv   
    73       COMPLEX(wp),DIMENSION(1) ::   y_fwfnow   
     73      COMPLEX(dp),DIMENSION(1) ::   y_fwfnow   
    7474      !!---------------------------------------------------------------------- 
    7575      ! 
    … …  
    180180            ! 
    181181!!gm   ===>>>>  lbc_lnk should be useless as all the computation is done over the whole domain ! 
    182             CALL lbc_lnk( 'sbcfwb', zerp_cor, 'T', 1. ) 
     182            CALL lbc_lnk( 'sbcfwb', zerp_cor, 'T', 1.0_wp ) 
    183183            ! 
    184184            emp(:,:) = emp(:,:) + zerp_cor(:,:) 
    … …  
    186186            erp(:,:) = erp(:,:) + zerp_cor(:,:) 
    187187            ! 
    188             IF( nprint == 1 .AND. lwp ) THEN                   ! control print 
     188            IF( lwp ) THEN                   ! control print 
    189189               IF( z_fwf < 0._wp ) THEN 
    190190                  WRITE(numout,*)'   z_fwf < 0' 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcice_cice.F90

    r12511 r13540  
    1212   USE oce             ! ocean dynamics and tracers 
    1313   USE dom_oce         ! ocean space and time domain 
     14# if ! defined key_qco 
    1415   USE domvvl 
     16# else 
     17   USE domqco 
     18# endif 
    1519   USE phycst, only : rcp, rho0, r1_rho0, rhos, rhoi 
    1620   USE in_out_manager  ! I/O manager 
    … …  
    213217! T point to U point 
    214218! T point to V point 
    215       DO_2D_10_10 
     219      DO_2D( 1, 0, 1, 0 ) 
    216220         fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1) 
    217221         fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1) 
    218222      END_2D 
    219223 
    220       CALL lbc_lnk_multi( 'sbcice_cice', fr_iu , 'U', 1.,  fr_iv , 'V', 1. ) 
     224      CALL lbc_lnk_multi( 'sbcice_cice', fr_iu , 'U', 1.0_wp,  fr_iv , 'V', 1.0_wp ) 
    221225 
    222226      ! set the snow+ice mass 
    … …  
    233237!!gm This should be put elsewhere....   (same remark for limsbc) 
    234238!!gm especially here it is assumed zstar coordinate, but it can be ztilde.... 
     239#if defined key_qco 
     240            IF( .NOT.ln_linssh )   CALL dom_qco_zgr( Kbb, Kmm, Kaa )   ! interpolation scale factor, depth and water column 
     241#else 
    235242            IF( .NOT.ln_linssh ) THEN 
    236243               ! 
    237244               DO jk = 1,jpkm1                     ! adjust initial vertical scale factors 
    238                   e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 
    239                   e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*tmask(:,:,1)/(ht_0(:,:) + 1.0 - tmask(:,:,1)) ) 
     245                  e3t(:,:,jk,Kmm) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kmm)*r1_ht_0(:,:)*tmask(:,:,jk) ) 
     246                  e3t(:,:,jk,Kbb) = e3t_0(:,:,jk)*( 1._wp + ssh(:,:,Kbb)*r1_ht_0(:,:)*tmask(:,:,jk) ) 
    240247               ENDDO 
    241248               e3t(:,:,:,Krhs) = e3t(:,:,:,Kbb) 
    … …  
    267274               END DO 
    268275            ENDIF 
     276#endif 
    269277         ENDIF 
    270278      ENDIF 
    … …  
    304312! x comp of wind stress (CI_1) 
    305313! U point to F point 
    306          DO_2D_10_11 
     314         DO_2D( 1, 0, 1, 1 ) 
    307315            ztmp(ji,jj) = 0.5 * (  fr_iu(ji,jj) * utau(ji,jj)      & 
    308316                                 + fr_iu(ji,jj+1) * utau(ji,jj+1) ) * fmask(ji,jj,1) 
    … …  
    312320! y comp of wind stress (CI_2) 
    313321! V point to F point 
    314          DO_2D_11_10 
     322         DO_2D( 1, 1, 1, 0 ) 
    315323            ztmp(ji,jj) = 0.5 * (  fr_iv(ji,jj) * vtau(ji,jj)      & 
    316324                                 + fr_iv(ji+1,jj) * vtau(ji+1,jj) ) * fmask(ji,jj,1) 
    … …  
    327335            qla_ice(:,:,1)= - ( emp_ice(:,:)+sprecip(:,:) ) * rLsub 
    328336! End of temporary code 
    329             DO_2D_11_11 
     337            DO_2D( 1, 1, 1, 1 ) 
    330338               IF(fr_i(ji,jj).eq.0.0) THEN 
    331339                  DO jl=1,ncat 
    … …  
    429437! x comp and y comp of surface ocean current 
    430438! U point to F point 
    431       DO_2D_10_11 
     439      DO_2D( 1, 0, 1, 1 ) 
    432440         ztmp(ji,jj)=0.5*(ssu_m(ji,jj)+ssu_m(ji,jj+1))*fmask(ji,jj,1) 
    433441      END_2D 
    … …  
    435443 
    436444! V point to F point 
    437       DO_2D_11_10 
     445      DO_2D( 1, 1, 1, 0 ) 
    438446         ztmp(ji,jj)=0.5*(ssv_m(ji,jj)+ssv_m(ji+1,jj))*fmask(ji,jj,1) 
    439447      END_2D 
    … …  
    459467! x comp and y comp of sea surface slope (on F points) 
    460468! T point to F point 
    461       DO_2D_10_10 
     469      DO_2D( 1, 0, 1, 0 ) 
    462470         ztmp(ji,jj)=0.5 * (  (zpice(ji+1,jj  )-zpice(ji,jj  )) * r1_e1u(ji,jj  )    & 
    463471            &               + (zpice(ji+1,jj+1)-zpice(ji,jj+1)) * r1_e1u(ji,jj+1)  ) * fmask(ji,jj,1) 
    … …  
    466474 
    467475! T point to F point 
    468       DO_2D_10_10 
     476      DO_2D( 1, 0, 1, 0 ) 
    469477         ztmp(ji,jj)=0.5 * (  (zpice(ji  ,jj+1)-zpice(ji  ,jj)) * r1_e2v(ji  ,jj)    & 
    470478            &               + (zpice(ji+1,jj+1)-zpice(ji+1,jj)) * r1_e2v(ji+1,jj)  ) *  fmask(ji,jj,1) 
    … …  
    495503      ss_iou(:,:)=0.0 
    496504! F point to U point 
    497       DO_2D_00_00 
     505      DO_2D( 0, 0, 0, 0 ) 
    498506         ss_iou(ji,jj) = 0.5 * ( ztmp1(ji,jj-1) + ztmp1(ji,jj) ) * umask(ji,jj,1) 
    499507      END_2D 
    500       CALL lbc_lnk( 'sbcice_cice', ss_iou , 'U', -1. ) 
     508      CALL lbc_lnk( 'sbcice_cice', ss_iou , 'U', -1.0_wp ) 
    501509 
    502510! y comp of ocean-ice stress  
    … …  
    505513! F point to V point 
    506514 
    507       DO_2D_10_00 
     515      DO_2D( 1, 0, 0, 0 ) 
    508516         ss_iov(ji,jj) = 0.5 * ( ztmp1(ji-1,jj) + ztmp1(ji,jj) ) * vmask(ji,jj,1) 
    509517      END_2D 
    510       CALL lbc_lnk( 'sbcice_cice', ss_iov , 'V', -1. ) 
     518      CALL lbc_lnk( 'sbcice_cice', ss_iov , 'V', -1.0_wp ) 
    511519 
    512520! x and y comps of surface stress 
    … …  
    561569      fmmflx(:,:) = ztmp1(:,:) !!Joakim edit 
    562570       
    563       CALL lbc_lnk_multi( 'sbcice_cice', emp , 'T', 1., sfx , 'T', 1. ) 
     571      CALL lbc_lnk_multi( 'sbcice_cice', emp , 'T', 1.0_wp, sfx , 'T', 1.0_wp ) 
    564572 
    565573! Solar penetrative radiation and non solar surface heat flux 
    … …  
    587595#endif 
    588596      qsr(:,:)=qsr(:,:)+ztmp1(:,:) 
    589       CALL lbc_lnk( 'sbcice_cice', qsr , 'T', 1. ) 
    590  
    591       DO_2D_11_11 
     597      CALL lbc_lnk( 'sbcice_cice', qsr , 'T', 1.0_wp ) 
     598 
     599      DO_2D( 1, 1, 1, 1 ) 
    592600         nfrzmlt(ji,jj)=MAX(nfrzmlt(ji,jj),0.0) 
    593601      END_2D 
    … …  
    600608      qns(:,:)=qns(:,:)+nfrzmlt(:,:)+ztmp1(:,:) 
    601609 
    602       CALL lbc_lnk( 'sbcice_cice', qns , 'T', 1. ) 
     610      CALL lbc_lnk( 'sbcice_cice', qns , 'T', 1.0_wp ) 
    603611 
    604612! Prepare for the following CICE time-step 
    … …  
    613621! T point to U point 
    614622! T point to V point 
    615       DO_2D_10_10 
     623      DO_2D( 1, 0, 1, 0 ) 
    616624         fr_iu(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji+1,jj))*umask(ji,jj,1) 
    617625         fr_iv(ji,jj)=0.5*(fr_i(ji,jj)+fr_i(ji,jj+1))*vmask(ji,jj,1) 
    618626      END_2D 
    619627 
    620       CALL lbc_lnk_multi( 'sbcice_cice', fr_iu , 'U', 1., fr_iv , 'V', 1. ) 
     628      CALL lbc_lnk_multi( 'sbcice_cice', fr_iu , 'U', 1.0_wp, fr_iv , 'V', 1.0_wp ) 
    621629 
    622630      ! set the snow+ice mass 
    … …  
    872880!        pcg(:,:)=0.0 
    873881         DO jn=1,jpnij 
    874             DO jj=nldjt(jn),nlejt(jn) 
    875                DO ji=nldit(jn),nleit(jn) 
     882            DO jj=njs0all(jn),nje0all(jn) 
     883               DO ji=nis0all(jn),nie0all(jn) 
    876884                  png2(ji+nimppt(jn)-1,jj+njmppt(jn)-1)=png(ji,jj,jn) 
    877885               ENDDO 
    … …  
    973981 
    974982      pn(:,:)=0.0 
    975       DO_2D_10_10 
     983      DO_2D( 1, 0, 1, 0 ) 
    976984         pn(ji,jj)=pc(ji+1-ji_off,jj+1-jj_off,1) 
    977985      END_2D 
    … …  
    9931001         png(:,:,:)=0.0 
    9941002         DO jn=1,jpnij 
    995             DO jj=nldjt(jn),nlejt(jn) 
    996                DO ji=nldit(jn),nleit(jn) 
     1003            DO jj=njs0all(jn),nje0all(jn) 
     1004               DO ji=nis0all(jn),nie0all(jn) 
    9971005                  png(ji,jj,jn)=pcg(ji+nimppt(jn)-1-ji_off,jj+njmppt(jn)-1-jj_off) 
    9981006               ENDDO 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcice_if.F90

    r12377 r13540  
    109109 
    110110         ! Flux and ice fraction computation 
    111          DO_2D_11_11 
     111         DO_2D( 1, 1, 1, 1 ) 
    112112            ! 
    113113            zt_fzp  = fr_i(ji,jj)                        ! freezing point temperature 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcmod.F90

    r12511 r13540  
    9999         &             nn_ice   , ln_ice_embd,                                       & 
    100100         &             ln_traqsr, ln_dm2dc ,                                         & 
    101          &             ln_rnf   , nn_fwb     , ln_ssr   , ln_apr_dyn,              & 
    102          &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc  , ln_stcor  ,     & 
     101         &             ln_rnf   , nn_fwb   , ln_ssr   , ln_apr_dyn,                  & 
     102         &             ln_wave  , ln_cdgw  , ln_sdw   , ln_tauwoc , ln_stcor  ,      & 
    103103         &             ln_tauw  , nn_lsm, nn_sdrift 
    104104      !!---------------------------------------------------------------------- 
    … …  
    120120      ncom_fsbc = nn_fsbc    ! make nn_fsbc available for lib_mpp 
    121121#endif 
    122       !                             !* overwrite namelist parameter using CPP key information 
    123 #if defined key_agrif 
    124       IF( Agrif_Root() ) THEN                ! AGRIF zoom (cf r1242: possibility to run without ice in fine grid) 
    125          IF( lk_si3  )   nn_ice      = 2 
    126          IF( lk_cice )   nn_ice      = 3 
    127       ENDIF 
    128 !!GS: TBD 
    129 !#else 
    130 !      IF( lk_si3  )   nn_ice      = 2 
    131 !      IF( lk_cice )   nn_ice      = 3 
     122#if ! defined key_si3 
     123      IF( nn_ice == 2 )    nn_ice = 0  ! without key key_si3 you cannot use si3... 
    132124#endif 
     125      ! 
    133126      ! 
    134127      IF(lwp) THEN                  !* Control print 
    … …  
    253246      ENDIF 
    254247      ! 
    255  
    256248      IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero 
    257249         IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case 
    … …  
    471463         ! A lbc_lnk is therefore needed to ensure reproducibility and restartability. 
    472464         ! see ticket #2113 for discussion about this lbc_lnk. 
    473          IF( .NOT. ln_passive_mode ) CALL lbc_lnk( 'sbcmod', emp, 'T', 1. ) ! ensure restartability with icebergs 
     465         IF( .NOT. ln_passive_mode ) CALL lbc_lnk( 'sbcmod', emp, 'T', 1.0_wp ) ! ensure restartability with icebergs 
    474466      ENDIF 
    475467 
    … …  
    486478!!$!RBbug do not understand why see ticket 667 
    487479!!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why. 
    488 !!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1. ) 
     480!!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1.0_wp ) 
    489481      IF( ll_wd ) THEN     ! If near WAD point limit the flux for now 
    490482         zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999  
    … …  
    517509            & iom_varid( numror, 'utau_b', ldstop = .FALSE. ) > 0 ) THEN 
    518510            IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields red in the restart file' 
    519             CALL iom_get( numror, jpdom_autoglo, 'utau_b', utau_b, ldxios = lrxios )   ! before i-stress  (U-point) 
    520             CALL iom_get( numror, jpdom_autoglo, 'vtau_b', vtau_b, ldxios = lrxios )   ! before j-stress  (V-point) 
    521             CALL iom_get( numror, jpdom_autoglo,  'qns_b',  qns_b, ldxios = lrxios )   ! before non solar heat flux (T-point) 
     511            CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b, ldxios = lrxios, cd_type = 'U', psgn = -1._wp )   ! before i-stress  (U-point) 
     512            CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b, ldxios = lrxios, cd_type = 'V', psgn = -1._wp )   ! before j-stress  (V-point) 
     513            CALL iom_get( numror, jpdom_auto,  'qns_b',  qns_b, ldxios = lrxios )   ! before non solar heat flux (T-point) 
    522514            ! The 3D heat content due to qsr forcing is treated in traqsr 
    523             ! CALL iom_get( numror, jpdom_autoglo, 'qsr_b' , qsr_b, ldxios = lrxios  ) ! before     solar heat flux (T-point) 
    524             CALL iom_get( numror, jpdom_autoglo, 'emp_b', emp_b, ldxios = lrxios  )    ! before     freshwater flux (T-point) 
     515            ! CALL iom_get( numror, jpdom_auto, 'qsr_b' , qsr_b, ldxios = lrxios  ) ! before     solar heat flux (T-point) 
     516            CALL iom_get( numror, jpdom_auto, 'emp_b', emp_b, ldxios = lrxios  )    ! before     freshwater flux (T-point) 
    525517            ! To ensure restart capability with 3.3x/3.4 restart files    !! to be removed in v3.6 
    526518            IF( iom_varid( numror, 'sfx_b', ldstop = .FALSE. ) > 0 ) THEN 
    527                CALL iom_get( numror, jpdom_autoglo, 'sfx_b', sfx_b, ldxios = lrxios )  ! before salt flux (T-point) 
     519               CALL iom_get( numror, jpdom_auto, 'sfx_b', sfx_b, ldxios = lrxios )  ! before salt flux (T-point) 
    528520            ELSE 
    529521               sfx_b (:,:) = sfx(:,:) 
    … …  
    573565      ENDIF 
    574566      ! 
    575       CALL iom_put( "utau", utau )   ! i-wind stress   (stress can be updated at each time step in sea-ice) 
    576       CALL iom_put( "vtau", vtau )   ! j-wind stress 
    577       ! 
    578567      IF(sn_cfctl%l_prtctl) THEN     ! print mean trends (used for debugging) 
    579568         CALL prt_ctl(tab2d_1=fr_i             , clinfo1=' fr_i     - : ', mask1=tmask ) 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcrnf.F90

    r12511 r13540  
    7272   !! * Substitutions 
    7373#  include "do_loop_substitute.h90" 
     74#  include "domzgr_substitute.h90" 
    7475   !!---------------------------------------------------------------------- 
    7576   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    159160            & iom_varid( numror, 'rnf_b', ldstop = .FALSE. ) > 0 ) THEN 
    160161            IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields red in the restart file', lrxios 
    161             CALL iom_get( numror, jpdom_autoglo, 'rnf_b', rnf_b, ldxios = lrxios )     ! before runoff 
    162             CALL iom_get( numror, jpdom_autoglo, 'rnf_hc_b', rnf_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before heat content of runoff 
    163             CALL iom_get( numror, jpdom_autoglo, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salinity content of runoff 
     162            CALL iom_get( numror, jpdom_auto, 'rnf_b', rnf_b, ldxios = lrxios )     ! before runoff 
     163            CALL iom_get( numror, jpdom_auto, 'rnf_hc_b', rnf_tsc_b(:,:,jp_tem), ldxios = lrxios )   ! before heat content of runoff 
     164            CALL iom_get( numror, jpdom_auto, 'rnf_sc_b', rnf_tsc_b(:,:,jp_sal), ldxios = lrxios )   ! before salinity content of runoff 
    164165         ELSE                                                   !* no restart: set from nit000 values 
    165166            IF(lwp) WRITE(numout,*) '          nit000-1 runoff forcing fields set to nit000' 
    … …  
    208209      IF( ln_rnf_depth .OR. ln_rnf_depth_ini ) THEN      !==   runoff distributed over several levels   ==! 
    209210         IF( ln_linssh ) THEN    !* constant volume case : just apply the runoff input flow 
    210             DO_2D_11_11 
     211            DO_2D( 1, 1, 1, 1 ) 
    211212               DO jk = 1, nk_rnf(ji,jj) 
    212213                  phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact * r1_rho0 / h_rnf(ji,jj) 
    … …  
    214215            END_2D 
    215216         ELSE                    !* variable volume case 
    216             DO_2D_11_11 
     217            DO_2D( 1, 1, 1, 1 )              ! update the depth over which runoffs are distributed 
    217218               h_rnf(ji,jj) = 0._wp 
    218                DO jk = 1, nk_rnf(ji,jj)                           ! recalculates h_rnf to be the depth in metres 
     219               DO jk = 1, nk_rnf(ji,jj)                             ! recalculates h_rnf to be the depth in metres 
    219220                  h_rnf(ji,jj) = h_rnf(ji,jj) + e3t(ji,jj,jk,Kmm)   ! to the bottom of the relevant grid box 
    220221               END DO 
    … …  
    353354         rn_dep_file = TRIM( cn_dir )//TRIM( sn_dep_rnf%clname ) 
    354355         IF( .NOT. sn_dep_rnf%ln_clim ) THEN   ;   WRITE(rn_dep_file, '(a,"_y",i4)' ) TRIM( rn_dep_file ), nyear    ! add year  
    355             IF( sn_dep_rnf%cltype == 'monthly' )   WRITE(rn_dep_file, '(a,"m",i2)'  ) TRIM( rn_dep_file ), nmonth   ! add month  
    356          ENDIF 
    357          CALL iom_open ( rn_dep_file, inum )                           ! open file 
    358          CALL iom_get  ( inum, jpdom_data, sn_dep_rnf%clvar, h_rnf )   ! read the river mouth array 
    359          CALL iom_close( inum )                                        ! close file 
     356            IF( sn_dep_rnf%clftyp == 'monthly' )   WRITE(rn_dep_file, '(a,"m",i2)'  ) TRIM( rn_dep_file ), nmonth   ! add month  
     357         ENDIF 
     358         CALL iom_open ( rn_dep_file, inum )                             ! open file 
     359         CALL iom_get  ( inum, jpdom_global, sn_dep_rnf%clvar, h_rnf )   ! read the river mouth array 
     360         CALL iom_close( inum )                                          ! close file 
    360361         ! 
    361362         nk_rnf(:,:) = 0                               ! set the number of level over which river runoffs are applied 
    362          DO_2D_11_11 
     363         DO_2D( 1, 1, 1, 1 ) 
    363364            IF( h_rnf(ji,jj) > 0._wp ) THEN 
    364365               jk = 2 
    … …  
    373374            ENDIF 
    374375         END_2D 
    375          DO_2D_11_11 
     376         DO_2D( 1, 1, 1, 1 )                           ! set the associated depth 
    376377            h_rnf(ji,jj) = 0._wp 
    377378            DO jk = 1, nk_rnf(ji,jj) 
    … …  
    390391         CALL iom_open( TRIM( sn_rnf%clname ), inum )    !  open runoff file 
    391392         nbrec = iom_getszuld( inum ) 
    392          zrnfcl(:,:,1) = 0._wp                                                          ! init the max to 0. in 1 
     393         zrnfcl(:,:,1) = 0._wp                                                            ! init the max to 0. in 1 
    393394         DO jm = 1, nbrec 
    394             CALL iom_get( inum, jpdom_data, TRIM( sn_rnf%clvar ), zrnfcl(:,:,2), jm )   ! read the value in 2 
    395             zrnfcl(:,:,1) = MAXVAL( zrnfcl(:,:,:), DIM=3 )                              ! store the maximum value in time in 1 
     395            CALL iom_get( inum, jpdom_global, TRIM( sn_rnf%clvar ), zrnfcl(:,:,2), jm )   ! read the value in 2 
     396            zrnfcl(:,:,1) = MAXVAL( zrnfcl(:,:,:), DIM=3 )                                ! store the maximum value in time in 1 
    396397         END DO 
    397398         CALL iom_close( inum ) 
    … …  
    403404         WHERE( zrnfcl(:,:,1) > 0._wp )  h_rnf(:,:) = zacoef * zrnfcl(:,:,1)   ! compute depth for all runoffs 
    404405         ! 
    405          DO_2D_11_11 
     406         DO_2D( 1, 1, 1, 1 )                ! take in account min depth of ocean rn_hmin 
    406407            IF( zrnfcl(ji,jj,1) > 0._wp ) THEN 
    407408               jk = mbkt(ji,jj) 
    … …  
    411412         ! 
    412413         nk_rnf(:,:) = 0                       ! number of levels on which runoffs are distributed 
    413          DO_2D_11_11 
     414         DO_2D( 1, 1, 1, 1 ) 
    414415            IF( zrnfcl(ji,jj,1) > 0._wp ) THEN 
    415416               jk = 2 
    … …  
    422423         END_2D 
    423424         ! 
    424          DO_2D_11_11 
     425         DO_2D( 1, 1, 1, 1 )                          ! set the associated depth 
    425426            h_rnf(ji,jj) = 0._wp 
    426427            DO jk = 1, nk_rnf(ji,jj) 
    … …  
    518519      cl_rnfile = TRIM( cn_dir )//TRIM( sn_cnf%clname ) 
    519520      IF( .NOT. sn_cnf%ln_clim ) THEN   ;   WRITE(cl_rnfile, '(a,"_y",i4)' ) TRIM( cl_rnfile ), nyear    ! add year 
    520          IF( sn_cnf%cltype == 'monthly' )   WRITE(cl_rnfile, '(a,"m",i2)'  ) TRIM( cl_rnfile ), nmonth   ! add month 
     521         IF( sn_cnf%clftyp == 'monthly' )   WRITE(cl_rnfile, '(a,"m",i2)'  ) TRIM( cl_rnfile ), nmonth   ! add month 
    521522      ENDIF 
    522523      ! 
    523524      ! horizontal mask (read in NetCDF file) 
    524       CALL iom_open ( cl_rnfile, inum )                           ! open file 
    525       CALL iom_get  ( inum, jpdom_data, sn_cnf%clvar, rnfmsk )    ! read the river mouth array 
    526       CALL iom_close( inum )                                      ! close file 
     525      CALL iom_open ( cl_rnfile, inum )                             ! open file 
     526      CALL iom_get  ( inum, jpdom_global, sn_cnf%clvar, rnfmsk )    ! read the river mouth array 
     527      CALL iom_close( inum )                                        ! close file 
    527528      ! 
    528529      IF( l_clo_rnf )   CALL clo_rnf( rnfmsk )   ! closed sea inflow set as river mouth 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcssm.F90

    r12377 r13540  
    3232   LOGICAL, SAVE ::   l_ssm_mean = .FALSE.   ! keep track of whether means have been read from restart file 
    3333    
     34#  include "domzgr_substitute.h90" 
    3435   !!---------------------------------------------------------------------- 
    3536   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    207208         IF( ln_rstart .AND. iom_varid( numror, 'nn_fsbc', ldstop = .FALSE. ) > 0 ) THEN 
    208209            l_ssm_mean = .TRUE. 
    209             CALL iom_get( numror               , 'nn_fsbc', zf_sbc, ldxios = lrxios )    ! sbc frequency of previous run 
    210             CALL iom_get( numror, jpdom_autoglo, 'ssu_m'  , ssu_m, ldxios = lrxios )    ! sea surface mean velocity    (U-point) 
    211             CALL iom_get( numror, jpdom_autoglo, 'ssv_m'  , ssv_m, ldxios = lrxios )    !   "         "    velocity    (V-point) 
    212             CALL iom_get( numror, jpdom_autoglo, 'sst_m'  , sst_m, ldxios = lrxios )    !   "         "    temperature (T-point) 
    213             CALL iom_get( numror, jpdom_autoglo, 'sss_m'  , sss_m, ldxios = lrxios )    !   "         "    salinity    (T-point) 
    214             CALL iom_get( numror, jpdom_autoglo, 'ssh_m'  , ssh_m, ldxios = lrxios )    !   "         "    height      (T-point) 
    215             CALL iom_get( numror, jpdom_autoglo, 'e3t_m'  , e3t_m, ldxios = lrxios )    ! 1st level thickness          (T-point) 
     210            CALL iom_get( numror            , 'nn_fsbc', zf_sbc,ldxios = lrxios )     ! sbc frequency of previous run 
     211            CALL iom_get( numror, jpdom_auto, 'ssu_m'  , ssu_m, ldxios = lrxios, cd_type = 'U', psgn = -1._wp )    ! sea surface mean velocity    (U-point) 
     212            CALL iom_get( numror, jpdom_auto, 'ssv_m'  , ssv_m, ldxios = lrxios, cd_type = 'V', psgn = -1._wp )    !   "         "    velocity    (V-point) 
     213            CALL iom_get( numror, jpdom_auto, 'sst_m'  , sst_m, ldxios = lrxios )    !   "         "    temperature (T-point) 
     214            CALL iom_get( numror, jpdom_auto, 'sss_m'  , sss_m, ldxios = lrxios )    !   "         "    salinity    (T-point) 
     215            CALL iom_get( numror, jpdom_auto, 'ssh_m'  , ssh_m, ldxios = lrxios )    !   "         "    height      (T-point) 
     216            CALL iom_get( numror, jpdom_auto, 'e3t_m'  , e3t_m, ldxios = lrxios )    ! 1st level thickness          (T-point) 
    216217            ! fraction of solar net radiation absorbed in 1st T level 
    217218            IF( iom_varid( numror, 'frq_m', ldstop = .FALSE. ) > 0 ) THEN 
    218                CALL iom_get( numror, jpdom_autoglo, 'frq_m'  , frq_m, ldxios = lrxios  ) 
     219               CALL iom_get( numror, jpdom_auto, 'frq_m'  , frq_m, ldxios = lrxios  ) 
    219220            ELSE 
    220221               frq_m(:,:) = 1._wp   ! default definition 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcssr.F90

    r12377 r13540  
    9595            ! 
    9696            IF( nn_sstr == 1 ) THEN                                   !* Temperature restoring term 
    97                DO_2D_11_11 
     97               DO_2D( 1, 1, 1, 1 ) 
    9898                  zqrp = rn_dqdt * ( sst_m(ji,jj) - sf_sst(1)%fnow(ji,jj,1) ) * tmask(ji,jj,1) 
    9999                  qns(ji,jj) = qns(ji,jj) + zqrp 
    … …  
    105105              ! use fraction of ice ( fr_i ) to adjust relaxation under ice if nn_sssr_ice .ne. 1 
    106106              ! n.b. coefice is initialised and fixed to 1._wp if nn_sssr_ice = 1 
    107                DO_2D_11_11 
     107               DO_2D( 1, 1, 1, 1 ) 
    108108                  SELECT CASE ( nn_sssr_ice ) 
    109109                    CASE ( 0 )    ;  coefice(ji,jj) = 1._wp - fr_i(ji,jj)              ! no/reduced damping under ice 
    … …  
    115115            IF( nn_sssr == 1 ) THEN                                   !* Salinity damping term (salt flux only (sfx)) 
    116116               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s] 
    117                DO_2D_11_11 
     117               DO_2D( 1, 1, 1, 1 ) 
    118118                  zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    119119                     &        *   coefice(ji,jj)            &      ! Optional control of damping under sea-ice 
    … …  
    126126               zsrp = rn_deds / rday                                  ! from [mm/day] to [kg/m2/s] 
    127127               zerp_bnd = rn_sssr_bnd / rday                          !       -              -     
    128                DO_2D_11_11 
     128               DO_2D( 1, 1, 1, 1 ) 
    129129                  zerp = zsrp * ( 1. - 2.*rnfmsk(ji,jj) )   &      ! No damping in vicinity of river mouths 
    130130                     &        *   coefice(ji,jj)            &      ! Optional control of damping under sea-ice 
    131131                     &        * ( sss_m(ji,jj) - sf_sss(1)%fnow(ji,jj,1) )   & 
    132132                     &        / MAX(  sss_m(ji,jj), 1.e-20   ) * tmask(ji,jj,1) 
    133                   IF( ln_sssr_bnd )   zerp = SIGN( 1., zerp ) * MIN( zerp_bnd, ABS(zerp) ) 
     133                  IF( ln_sssr_bnd )   zerp = SIGN( 1.0_wp, zerp ) * MIN( zerp_bnd, ABS(zerp) ) 
    134134                  emp(ji,jj) = emp (ji,jj) + zerp 
    135135                  qns(ji,jj) = qns(ji,jj) - zerp * rcp * sst_m(ji,jj) 

  • NEMO/branches/2020/r12377_ticket2386/src/OCE/SBC/sbcwave.F90

    r12377 r13540  
    7373   !! * Substitutions 
    7474#  include "do_loop_substitute.h90" 
     75#  include "domzgr_substitute.h90" 
    7576   !!---------------------------------------------------------------------- 
    7677   !! NEMO/OCE 4.0 , NEMO Consortium (2018) 
    … …  
    112113      IF( ll_st_bv_li ) THEN   ! (Eq. (19) in Breivik et al. (2014) ) 
    113114         zfac = 2.0_wp * rpi / 16.0_wp 
    114          DO_2D_11_11 
     115         DO_2D( 1, 1, 1, 1 ) 
    115116            ! Stokes drift velocity estimated from Hs and Tmean 
    116117            ztransp = zfac * hsw(ji,jj)*hsw(ji,jj) / MAX( wmp(ji,jj), 0.0000001_wp ) 
    … …  
    120121            zk_t(ji,jj) = ABS( tsd2d(ji,jj) ) / MAX( ABS( 5.97_wp*ztransp ), 0.0000001_wp ) 
    121122         END_2D 
    122          DO_2D_10_10 
     123         DO_2D( 1, 0, 1, 0 )          ! exp. wave number & Stokes drift velocity at u- & v-points 
    123124            zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 
    124125            zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 
    … …  
    128129         END_2D 
    129130      ELSE IF( ll_st_peakfr ) THEN    ! peak wave number calculated from the peak frequency received by the wave model 
    130          DO_2D_11_11 
     131         DO_2D( 1, 1, 1, 1 ) 
    131132            zk_t(ji,jj) = ( 2.0_wp * rpi * wfreq(ji,jj) ) * ( 2.0_wp * rpi * wfreq(ji,jj) ) / grav 
    132133         END_2D 
    133          DO_2D_10_10 
     134         DO_2D( 1, 0, 1, 0 ) 
    134135            zk_u(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji+1,jj) ) 
    135136            zk_v(ji,jj) = 0.5_wp * ( zk_t(ji,jj) + zk_t(ji,jj+1) ) 
    … …  
    142143      !                       !==  horizontal Stokes Drift 3D velocity  ==! 
    143144      IF( ll_st_bv2014 ) THEN 
    144          DO_3D_00_00( 1, jpkm1 ) 
     145         DO_3D( 0, 0, 0, 0, 1, jpkm1 ) 
    145146            zdep_u = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji+1,jj,jk,Kmm) ) 
    146147            zdep_v = 0.5_wp * ( gdept(ji,jj,jk,Kmm) + gdept(ji,jj+1,jk,Kmm) ) 
    … …  
    157158      ELSE IF( ll_st_li2017 .OR. ll_st_peakfr ) THEN 
    158159         ALLOCATE( zstokes_psi_u_top(jpi,jpj), zstokes_psi_v_top(jpi,jpj) ) 
    159          DO_2D_10_10 
     160         DO_2D( 1, 0, 1, 0 ) 
    160161            zstokes_psi_u_top(ji,jj) = 0._wp 
    161162            zstokes_psi_v_top(ji,jj) = 0._wp 
    … …  
    163164         zsqrtpi = SQRT(rpi) 
    164165         z_two_thirds = 2.0_wp / 3.0_wp 
    165          DO_3D_00_00( 1, jpkm1 ) 
     166         DO_3D( 0, 0, 0, 0, 1, jpkm1 )       ! exp. wave number & Stokes drift velocity at u- & v-points 
    166167            zbot_u = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji+1,jj,jk+1,Kmm) )  ! 2 * bottom depth 
    167168            zbot_v = ( gdepw(ji,jj,jk+1,Kmm) + gdepw(ji,jj+1,jk+1,Kmm) )  ! 2 * bottom depth 
    … …  
    198199      ENDIF 
    199200 
    200       CALL lbc_lnk_multi( 'sbcwave', usd, 'U', -1., vsd, 'V', -1. ) 
     201      CALL lbc_lnk_multi( 'sbcwave', usd, 'U', -1.0_wp, vsd, 'V', -1.0_wp ) 
    201202 
    202203      ! 
    203204      !                       !==  vertical Stokes Drift 3D velocity  ==! 
    204205      ! 
    205       DO_3D_01_01( 1, jpkm1 ) 
     206      DO_3D( 0, 1, 0, 1, 1, jpkm1 )    ! Horizontal e3*divergence 
    206207         ze3divh(ji,jj,jk) = (  e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * usd(ji  ,jj,jk)    & 
    207208            &                 - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * usd(ji-1,jj,jk)    & 
    208209            &                 + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * vsd(ji,jj  ,jk)    & 
    209             &                 - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vsd(ji,jj-1,jk)  ) * r1_e1e2t(ji,jj) 
     210            &                 - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * vsd(ji,jj-1,jk)  ) & 
     211            &                * r1_e1e2t(ji,jj) 
    210212      END_3D 
    211213      ! 
    212 #if defined key_agrif 
    213       IF( .NOT. Agrif_Root() ) THEN 
    214          IF( nbondi == -1 .OR. nbondi == 2 )   ze3divh( 2:nbghostcells+1,:        ,:) = 0._wp      ! west 
    215          IF( nbondi ==  1 .OR. nbondi == 2 )   ze3divh( nlci-nbghostcells:nlci-1,:,:) = 0._wp      ! east 
    216          IF( nbondj == -1 .OR. nbondj == 2 )   ze3divh( :,2:nbghostcells+1        ,:) = 0._wp      ! south 
    217          IF( nbondj ==  1 .OR. nbondj == 2 )   ze3divh( :,nlcj-nbghostcells:nlcj-1,:) = 0._wp      ! north 
    218       ENDIF 
    219 #endif 
    220       ! 
    221       CALL lbc_lnk( 'sbcwave', ze3divh, 'T', 1. ) 
     214      CALL lbc_lnk( 'sbcwave', ze3divh, 'T', 1.0_wp ) 
    222215      ! 
    223216      IF( ln_linssh ) THEN   ;   ik = 1   ! none zero velocity through the sea surface 
    … …  
    270263      ! 
    271264      IF( ln_tauw ) THEN 
    272          DO_2D_10_10 
     265         DO_2D( 1, 0, 1, 0 ) 
    273266            ! Stress components at u- & v-points 
    274267            utau(ji,jj) = 0.5_wp * ( tauw_x(ji,jj) + tauw_x(ji+1,jj) ) 
    … …  
    278271            taum(ji,jj) = SQRT( tauw_x(ji,jj)*tauw_x(ji,jj) + tauw_y(ji,jj)*tauw_y(ji,jj) ) 
    279272         END_2D 
    280          CALL lbc_lnk_multi( 'sbcwave', utau(:,:), 'U', -1. , vtau(:,:), 'V', -1. , taum(:,:) , 'T', -1. ) 
     273         CALL lbc_lnk_multi( 'sbcwave', utau(:,:), 'U', -1.0_wp , vtau(:,:), 'V', -1.0_wp , taum(:,:) , 'T', -1.0_wp ) 
    281274      ENDIF 
    282275      ! 
Note: See TracChangeset for help on using the changeset viewer.