Changeset 11963


Ignore:
Timestamp:
2019-11-26T12:08:01+01:00 (8 weeks ago)
Author:
laurent
Message:

More accurate comments/info, better syntax, simplifications, etc

Location:
NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC
Files:
20 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/cpl_oasis3.F90

    r10582 r11963  
    114114      !------------------------------------------------------------------ 
    115115      CALL oasis_init_comp ( ncomp_id, TRIM(cd_modname), nerror ) 
    116       IF ( nerror /= OASIS_Ok ) & 
     116      IF( nerror /= OASIS_Ok ) & 
    117117         CALL oasis_abort (ncomp_id, 'cpl_init', 'Failure in oasis_init_comp') 
    118118 
     
    122122 
    123123      CALL oasis_get_localcomm ( kl_comm, nerror ) 
    124       IF ( nerror /= OASIS_Ok ) & 
     124      IF( nerror /= OASIS_Ok ) & 
    125125         CALL oasis_abort (ncomp_id, 'cpl_init','Failure in oasis_get_localcomm' ) 
    126126      ! 
     
    149149 
    150150      ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 
    151       IF ( ltmp_wapatch ) THEN 
     151      IF( ltmp_wapatch ) THEN 
    152152         nldi_save = nldi   ;   nlei_save = nlei 
    153153         nldj_save = nldj   ;   nlej_save = nlej 
     
    217217      ! 
    218218      DO ji = 1, ksnd 
    219          IF ( ssnd(ji)%laction ) THEN 
     219         IF( ssnd(ji)%laction ) THEN 
    220220 
    221221            IF( ssnd(ji)%nct > nmaxcat ) THEN 
     
    228228               DO jm = 1, kcplmodel 
    229229 
    230                   IF ( ssnd(ji)%nct .GT. 1 ) THEN 
     230                  IF( ssnd(ji)%nct .GT. 1 ) THEN 
    231231                     WRITE(cli2,'(i2.2)') jc 
    232232                     zclname = TRIM(ssnd(ji)%clname)//'_cat'//cli2 
     
    234234                     zclname = ssnd(ji)%clname 
    235235                  ENDIF 
    236                   IF ( kcplmodel  > 1 ) THEN 
     236                  IF( kcplmodel  > 1 ) THEN 
    237237                     WRITE(cli2,'(i2.2)') jm 
    238238                     zclname = 'model'//cli2//'_'//TRIM(zclname) 
     
    241241                  IF( agrif_fixed() /= 0 ) THEN  
    242242                     zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname) 
    243                   END IF 
     243                  ENDIF 
    244244#endif 
    245245                  IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_Out 
    246246                  CALL oasis_def_var (ssnd(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 1 /),   & 
    247247                     &                OASIS_Out          , ishape , OASIS_REAL, nerror ) 
    248                   IF ( nerror /= OASIS_Ok ) THEN 
     248                  IF( nerror /= OASIS_Ok ) THEN 
    249249                     WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname) 
    250250                     CALL oasis_abort ( ssnd(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' ) 
     
    262262      ! 
    263263      DO ji = 1, krcv 
    264          IF ( srcv(ji)%laction ) THEN  
     264         IF( srcv(ji)%laction ) THEN  
    265265             
    266266            IF( srcv(ji)%nct > nmaxcat ) THEN 
     
    273273               DO jm = 1, kcplmodel 
    274274                   
    275                   IF ( srcv(ji)%nct .GT. 1 ) THEN 
     275                  IF( srcv(ji)%nct .GT. 1 ) THEN 
    276276                     WRITE(cli2,'(i2.2)') jc 
    277277                     zclname = TRIM(srcv(ji)%clname)//'_cat'//cli2 
     
    279279                     zclname = srcv(ji)%clname 
    280280                  ENDIF 
    281                   IF ( kcplmodel  > 1 ) THEN 
     281                  IF( kcplmodel  > 1 ) THEN 
    282282                     WRITE(cli2,'(i2.2)') jm 
    283283                     zclname = 'model'//cli2//'_'//TRIM(zclname) 
     
    286286                  IF( agrif_fixed() /= 0 ) THEN  
    287287                     zclname=TRIM(Agrif_CFixed())//'_'//TRIM(zclname) 
    288                   END IF 
     288                  ENDIF 
    289289#endif 
    290290                  IF( ln_ctl ) WRITE(numout,*) "Define", ji, jc, jm, " "//TRIM(zclname), " for ", OASIS_In 
    291291                  CALL oasis_def_var (srcv(ji)%nid(jc,jm), zclname, id_part   , (/ 2, 1 /),   & 
    292292                     &                OASIS_In           , ishape , OASIS_REAL, nerror ) 
    293                   IF ( nerror /= OASIS_Ok ) THEN 
     293                  IF( nerror /= OASIS_Ok ) THEN 
    294294                     WRITE(numout,*) 'Failed to define transient ', ji, jc, jm, " "//TRIM(zclname) 
    295295                     CALL oasis_abort ( srcv(ji)%nid(jc,jm), 'cpl_define', 'Failure in oasis_def_var' ) 
     
    310310      IF( nerror /= OASIS_Ok )   CALL oasis_abort ( ncomp_id, 'cpl_define', 'Failure in oasis_enddef') 
    311311      ! 
    312       IF ( ltmp_wapatch ) THEN 
     312      IF( ltmp_wapatch ) THEN 
    313313         nldi = nldi_save   ;   nlei = nlei_save 
    314314         nldj = nldj_save   ;   nlej = nlej_save 
     
    332332      !!-------------------------------------------------------------------- 
    333333      ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 
    334       IF ( ltmp_wapatch ) THEN 
     334      IF( ltmp_wapatch ) THEN 
    335335         nldi_save = nldi   ;   nlei_save = nlei 
    336336         nldj_save = nldj   ;   nlej_save = nlej 
     
    349349               CALL oasis_put ( ssnd(kid)%nid(jc,jm), kstep, pdata(nldi:nlei, nldj:nlej,jc), kinfo ) 
    350350                
    351                IF ( ln_ctl ) THEN         
    352                   IF ( kinfo == OASIS_Sent     .OR. kinfo == OASIS_ToRest .OR.   & 
     351               IF( ln_ctl ) THEN         
     352                  IF( kinfo == OASIS_Sent     .OR. kinfo == OASIS_ToRest .OR.   & 
    353353                     & kinfo == OASIS_SentOut  .OR. kinfo == OASIS_ToRestOut ) THEN 
    354354                     WRITE(numout,*) '****************' 
     
    368368         ENDDO 
    369369      ENDDO 
    370       IF ( ltmp_wapatch ) THEN 
     370      IF( ltmp_wapatch ) THEN 
    371371         nldi = nldi_save   ;   nlei = nlei_save 
    372372         nldj = nldj_save   ;   nlej = nlej_save 
     
    393393      !!-------------------------------------------------------------------- 
    394394      ! patch to restore wraparound rows in cpl_send, cpl_rcv, cpl_define 
    395       IF ( ltmp_wapatch ) THEN 
     395      IF( ltmp_wapatch ) THEN 
    396396         nldi_save = nldi   ;   nlei_save = nlei 
    397397         nldj_save = nldj   ;   nlej_save = nlej 
     
    403403      ! 
    404404      DO jc = 1, srcv(kid)%nct 
    405          IF ( ltmp_wapatch ) THEN 
     405         IF( ltmp_wapatch ) THEN 
    406406            IF( nimpp           ==      1 ) nldi = 1 
    407407            IF( nimpp + jpi - 1 == jpiglo ) nlei = jpi 
     
    420420                  &        kinfo == OASIS_RecvOut .OR. kinfo == OASIS_FromRestOut 
    421421                
    422                IF ( ln_ctl )   WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm) 
     422               IF( ln_ctl )   WRITE(numout,*) "llaction, kinfo, kstep, ivarid: " , llaction, kinfo, kstep, srcv(kid)%nid(jc,jm) 
    423423                
    424                IF ( llaction ) THEN 
     424               IF( llaction ) THEN 
    425425                   
    426426                  kinfo = OASIS_Rcv 
     
    432432                  ENDIF 
    433433                   
    434                   IF ( ln_ctl ) THEN         
     434                  IF( ln_ctl ) THEN         
    435435                     WRITE(numout,*) '****************' 
    436436                     WRITE(numout,*) 'oasis_get: Incoming ', srcv(kid)%clname 
     
    450450         ENDDO 
    451451 
    452          IF ( ltmp_wapatch ) THEN 
     452         IF( ltmp_wapatch ) THEN 
    453453            nldi = nldi_save   ;   nlei = nlei_save 
    454454            nldj = nldj_save   ;   nlej = nlej_save 
     
    483483      ! 
    484484      DO ji = 1, nsnd 
    485          IF (ssnd(ji)%laction ) THEN 
     485         IF(ssnd(ji)%laction ) THEN 
    486486            DO jm = 1, ncplmodel 
    487487               IF( ssnd(ji)%nid(1,jm) /= -1 ) THEN 
     
    495495      ENDDO 
    496496      DO ji = 1, nrcv 
    497          IF (srcv(ji)%laction ) THEN 
     497         IF(srcv(ji)%laction ) THEN 
    498498            DO jm = 1, ncplmodel 
    499499               IF( srcv(ji)%nid(1,jm) /= -1 ) THEN 
     
    529529      ! 
    530530      DEALLOCATE( exfld ) 
    531       IF (nstop == 0) THEN 
     531      IF(nstop == 0) THEN 
    532532         CALL oasis_terminate( nerror )          
    533533      ELSE 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/cyclone.F90

    r10068 r11963  
    137137            zhemi = SIGN( 1. , zrlat ) 
    138138            zinfl = 15.* rad                             ! clim inflow angle in Tropical Cyclones 
    139          IF ( vortex == 0 ) THEN 
     139         IF( vortex == 0 ) THEN 
    140140 
    141141            ! Vortex Holland reconstruct wind at each lon-lat position 
     
    157157                     &              + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) ) 
    158158 
    159                  IF (zdist < zrout2) THEN ! calculation of wind only to a given max radius 
     159                 IF(zdist < zrout2) THEN ! calculation of wind only to a given max radius 
    160160                  ! shape of the wind profile 
    161161                  zztmp = ( zrmw / ( zdist + 1.e-12 ) )**zb 
    162162                  zztmp =  zvmax * SQRT( zztmp * EXP(1. - zztmp) )     
    163163 
    164                   IF (zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 
     164                  IF(zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 
    165165                     zztmp = zztmp * ( (zrout2-zdist)*1.e-6 ) 
    166166                  ENDIF 
    167167 
    168168                  ! !!! KILL EQ WINDS 
    169                   ! IF (SIGN( 1. , zrlat ) /= zhemi) THEN 
     169                  ! IF(SIGN( 1. , zrlat ) /= zhemi) THEN 
    170170                  !    zztmp = 0.                              ! winds in other hemisphere 
    171                   !    IF (ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S 
    172                   ! ENDIF 
    173                   ! IF (ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 
     171                  !    IF(ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S 
     172                  ! ENDIF 
     173                  ! IF(ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 
    174174                  !    zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) )  
    175175                  !    !linear to zero between 10 and 5 
     
    177177                  ! !!! / KILL EQ 
    178178 
    179                   IF (ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 
     179                  IF(ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 
    180180 
    181181                  zwnd_t =   COS( zinfl ) * zztmp     
     
    196196            END DO 
    197197          
    198          ELSE IF ( vortex == 1 ) THEN 
     198         ELSE IF( vortex == 1 ) THEN 
    199199 
    200200            ! Vortex Willoughby reconstruct wind at each lon-lat position 
     
    206206            zn   =   2.1340 + 0.0077*zvmax - 0.4522*LOG(zrmw/1000.) - 0.0038*ABS( ztct(jtc,jp_lat) )             
    207207            zA   =   0.5913 + 0.0029*zvmax - 0.1361*LOG(zrmw/1000.) - 0.0042*ABS( ztct(jtc,jp_lat) )   
    208             IF (zA < 0) THEN  
     208            IF(zA < 0) THEN  
    209209               zA=0 
    210210            ENDIF            
     
    218218                     &              + COS( zrlat ) * COS( zzrgphi ) * COS( zzrglam ) ) 
    219219 
    220                  IF (zdist < zrout2) THEN ! calculation of wind only to a given max radius 
     220                 IF(zdist < zrout2) THEN ! calculation of wind only to a given max radius 
    221221                
    222222                  ! shape of the wind profile                      
    223                   IF (zdist <= zrmw) THEN     ! inside the Radius of Maximum Wind 
     223                  IF(zdist <= zrmw) THEN     ! inside the Radius of Maximum Wind 
    224224                     zztmp  = zvmax * (zdist/zrmw)**zn 
    225225                  ELSE  
     
    227227                  ENDIF 
    228228 
    229                   IF (zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 
     229                  IF(zdist > zrout1) THEN ! bring to zero between r_out1 and r_out2 
    230230                     zztmp = zztmp * ( (zrout2-zdist)*1.e-6 ) 
    231231                  ENDIF 
    232232 
    233233                  ! !!! KILL EQ WINDS 
    234                   ! IF (SIGN( 1. , zrlat ) /= zhemi) THEN 
     234                  ! IF(SIGN( 1. , zrlat ) /= zhemi) THEN 
    235235                  !    zztmp = 0.                              ! winds in other hemisphere 
    236                   !    IF (ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S 
    237                   ! ENDIF 
    238                   ! IF (ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 
     236                  !    IF(ABS(gphit(ji,jj)) <= 5.) zztmp=0.   ! kill between 5N-5S 
     237                  ! ENDIF 
     238                  ! IF(ABS(gphit(ji,jj)) <= 10. .and. ABS(gphit(ji,jj)) > 5.) THEN 
    239239                  !    zztmp = zztmp * ( 1./5. * (ABS(gphit(ji,jj)) - 5.) )  
    240240                  !    !linear to zero between 10 and 5 
     
    242242                  ! !!! / KILL EQ 
    243243 
    244                   IF (ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 
     244                  IF(ABS(gphit(ji,jj)) >= 55.) zztmp = 0. ! kill weak spurious winds at high latitude 
    245245 
    246246                  zwnd_t =   COS( zinfl ) * zztmp     
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/fldread.F90

    r11831 r11963  
    167167      IF( PRESENT(kit) )   ll_firstcall = ll_firstcall .and. kit == 1 
    168168 
    169       IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc 
     169      IF( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc 
    170170      ELSE                                      ;   it_offset = 0 
    171171      ENDIF 
     
    389389         ENDIF 
    390390         ! 
    391          IF ( sdjf%cltype(1:4) == 'week' ) THEN 
     391         IF( sdjf%cltype(1:4) == 'week' ) THEN 
    392392            isec_week = isec_week + ksec_week( sdjf%cltype(6:8) )   ! second since the beginning of the week 
    393393            llprevmth = isec_week > nsec_month                      ! longer time since the beginning of the week than the month 
     
    464464      ENDIF 
    465465      ! 
    466       IF ( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc 
     466      IF( nn_components == jp_iam_sas ) THEN   ;   it_offset = nn_fsbc 
    467467      ELSE                                      ;   it_offset = 0 
    468468      ENDIF 
     
    656656            ENDIF 
    657657         CASE DEFAULT 
    658             IF (lk_c1d .AND. lmoor ) THEN 
     658            IF(lk_c1d .AND. lmoor ) THEN 
    659659               IF( sdjf%ln_tint ) THEN 
    660660                  CALL iom_get( sdjf%num, jpdom_unknown, sdjf%clvar, sdjf%fdta(2,2,:,2), sdjf%nrec_a(1) ) 
     
    10711071         imonth = kmonth 
    10721072         iday = kday 
    1073          IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week 
     1073         IF( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week 
    10741074            isec_week = ksec_week( sdjf%cltype(6:8) )- (86400 * 8 )   
    10751075            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month 
     
    10801080         ENDIF 
    10811081      ELSE                                                  ! use current day values 
    1082          IF ( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week 
     1082         IF( sdjf%cltype(1:4) == 'week' ) THEN             ! find the day of the beginning of the week 
    10831083            isec_week  = ksec_week( sdjf%cltype(6:8) )      ! second since the beginning of the week 
    10841084            llprevmth  = isec_week > nsec_month             ! longer time since beginning of the week than the month 
     
    13181318 
    13191319      !! get dimensions 
    1320       IF ( SIZE(sd%fnow, 3) > 1 ) THEN 
     1320      IF( SIZE(sd%fnow, 3) > 1 ) THEN 
    13211321         ALLOCATE( ddims(4) ) 
    13221322      ELSE 
     
    13311331 
    13321332      CALL iom_open ( sd%wgtname, inum )   ! interpolation weights 
    1333       IF ( inum > 0 ) THEN 
     1333      IF( inum > 0 ) THEN 
    13341334 
    13351335         !! determine whether we have an east-west cyclic grid 
     
    16631663      END DO 
    16641664 
    1665       IF (ref_wgts(kw)%numwgt .EQ. 16) THEN 
     1665      IF(ref_wgts(kw)%numwgt .EQ. 16) THEN 
    16661666 
    16671667        !! fix up halo points that we couldnt read from file 
     
    17461746         END DO 
    17471747         ! 
    1748       END IF 
     1748      ENDIF 
    17491749      ! 
    17501750   END SUBROUTINE fld_interp 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcapr.F90

    r11831 r11963  
    103103      ! 
    104104      !                                            !* control check 
    105       IF ( ln_apr_obc  ) THEN 
     105      IF( ln_apr_obc  ) THEN 
    106106         IF(lwp) WRITE(numout,*) '         Inverse barometer added to OBC ssh data' 
    107107      ENDIF 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk.F90

    r11962 r11963  
    259259      CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_init', 'surface boundary condition -- bulk formulae', 'namsbc_blk' ) 
    260260      ! 
    261       IF ( ln_wave ) THEN 
     261      IF( ln_wave ) THEN 
    262262         !Activated wave module but neither drag nor stokes drift activated 
    263          IF ( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN 
     263         IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) )   THEN 
    264264            CALL ctl_stop( 'STOP',  'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' ) 
    265265            !drag coefficient read from wave model definable only with mfs bulk formulae and core 
    266          ELSEIF (ln_cdgw .AND. .NOT. ln_NCAR )       THEN 
     266         ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR )       THEN 
    267267            CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae') 
    268          ELSEIF (ln_stcor .AND. .NOT. ln_sdw)                             THEN 
     268         ELSEIF(ln_stcor .AND. .NOT. ln_sdw)                             THEN 
    269269            CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T') 
    270270         ENDIF 
    271271      ELSE 
    272          IF ( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                & 
     272         IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor )                & 
    273273            &   CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ',    & 
    274274            &                  'with drag coefficient (ln_cdgw =T) '  ,                        & 
     
    481481      zsq(:,:) = rdct_qsat_salt * q_sat( zst(:,:), sf(jp_slp)%fnow(:,:,1) ) 
    482482 
    483       IF ( ln_skin_cs .OR. ln_skin_wl ) THEN 
     483      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    484484         zqsb(:,:) = zst(:,:) !LB: using array zqsb to backup original zst before skin action 
    485485         zqla(:,:) = zsq(:,:) !LB: using array zqla to backup original zsq before skin action 
    486       END IF 
     486      ENDIF 
    487487 
    488488      !LB: 
     
    492492         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity! 
    493493      CASE( np_humi_dpt ) 
    494          !IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 
     494         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of d_air and slp !' !LBrm 
    495495         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    496496      CASE( np_humi_rlh ) 
    497          !IF (lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 
     497         !IF(lwp) WRITE(numout,*) ' *** blk_oce => computing q_air out of RH, t_air and slp !' !LBrm 
    498498         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file 
    499499      END SELECT 
     
    506506 
    507507 
    508       IF ( ln_skin_cs .OR. ln_skin_wl ) THEN 
     508      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    509509 
    510510         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
     
    543543         tsk(:,:) = zst(:,:) 
    544544 
    545       ELSE !IF ( ln_skin_cs .OR. ln_skin_wl ) 
     545      ELSE !IF( ln_skin_cs .OR. ln_skin_wl ) 
    546546 
    547547         SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point 
     
    567567         END SELECT 
    568568 
    569       END IF ! IF ( ln_skin_cs .OR. ln_skin_wl ) 
     569      ENDIF ! IF( ln_skin_cs .OR. ln_skin_wl ) 
    570570 
    571571      !!      CALL iom_put( "Cd_oce", Cd_atm)  ! output value of pure ocean-atm. transfer coef. 
     
    576576         t_zu(:,:) = ztpot(:,:) 
    577577         q_zu(:,:) = zqair(:,:) 
    578       END IF 
     578      ENDIF 
    579579 
    580580 
     
    660660#endif 
    661661      ! 
    662       !!#LB: NO WHY???? IF ( nn_ice == 0 ) THEN 
     662      !!#LB: NO WHY???? IF( nn_ice == 0 ) THEN 
    663663      CALL iom_put( "rho_air"  ,   rhoa )                 ! output air density (kg/m^3) !#LB 
    664664      CALL iom_put( "qlw_oce"  ,   zqlw )                 ! output downward longwave heat over the ocean 
     
    674674      CALL iom_put( 'snowpre', sprecip )                 ! Snow 
    675675      CALL iom_put( 'precip' , tprecip )                 ! Total precipitation 
    676       IF ( ln_skin_cs .OR. ln_skin_wl ) THEN 
     676      IF( ln_skin_cs .OR. ln_skin_wl ) THEN 
    677677         CALL iom_put( "t_skin" ,  (zst - rt0) * tmask(:,:,1) )           ! T_skin in Celsius 
    678678         CALL iom_put( "dt_skin" , (zst - pst - rt0) * tmask(:,:,1) )     ! T_skin - SST temperature difference... 
    679       END IF 
     679      ENDIF 
    680680      !!#LB. ENDIF 
    681681      ! 
     
    819819         zqair(:,:) =        sf(jp_humi)%fnow(:,:,1)      ! what we read in file is already a spec. humidity! 
    820820      CASE( np_humi_dpt ) 
    821          !IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm 
     821         !IF(lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of d_air and slp !' !LBrm 
    822822         zqair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) 
    823823      CASE( np_humi_rlh ) 
    824          !IF (lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm 
     824         !IF(lwp) WRITE(numout,*) ' *** blk_ice_flx => computing q_air out of RH, t_air and slp !' !LBrm 
    825825         zqair(:,:) = q_air_rh( 0.01_wp*sf(jp_humi)%fnow(:,:,1), sf(jp_tair)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) ) !LB: 0.01 => RH is % percent in file 
    826826      END SELECT 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p0.F90

    r11962 r11963  
    6868      INTEGER :: ierr 
    6969      !!--------------------------------------------------------------------- 
    70       IF ( l_use_wl ) THEN 
     70      IF( l_use_wl ) THEN 
    7171         ierr = 0 
    7272         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) 
     
    7676         dT_wl(:,:)  = 0._wp 
    7777         Hz_wl(:,:)  = Hwl_max 
    78       END IF 
    79       IF ( l_use_cs ) THEN 
     78      ENDIF 
     79      IF( l_use_cs ) THEN 
    8080         ierr = 0 
    8181         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 
    8282         IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_COARE3P0_INIT => allocation of dT_cs failed!' ) 
    8383         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction 
    84       END IF 
     84      ENDIF 
    8585   END SUBROUTINE sbcblk_algo_coare3p0_init 
    8686 
     
    190190      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p0@sbcblk_algo_coare3p0' 
    191191      !!---------------------------------------------------------------------------------- 
    192       IF ( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P0_INIT(l_use_cs, l_use_wl) 
     192      IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P0_INIT(l_use_cs, l_use_wl) 
    193193 
    194194      l_zt_equal_zu = .FALSE. 
     
    197197 
    198198      !! Initializations for cool skin and warm layer: 
    199       IF ( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
     199      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
    200200         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' ) 
    201201 
    202       IF ( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
     202      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
    203203         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' ) 
    204204 
    205       IF ( l_use_cs .OR. l_use_wl ) THEN 
     205      IF( l_use_cs .OR. l_use_wl ) THEN 
    206206         ALLOCATE ( zsst(jpi,jpj) ) 
    207207         zsst = T_s ! backing up the bulk SST 
    208208         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction 
    209209         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
    210       END IF 
     210      ENDIF 
    211211 
    212212 
     
    265265         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    266266         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    267       END IF 
     267      ENDIF 
    268268 
    269269      !! ITERATION BLOCK 
     
    288288            zeta_t = zt*ztmp0 
    289289            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t ) 
    290          END IF 
     290         ENDIF 
    291291 
    292292         !! Adjustment the wind at 10m (not needed in the current algo form): 
    293          !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u)) 
     293         !IF( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u)) 
    294294 
    295295         !! Roughness lengthes z0, z0t (z0q = z0t) : 
     
    299299 
    300300         ztmp1 = ( znu_a / (z0*u_star) )**0.6_wp    ! (1./Re_r)^0.72 (Re_r: roughness Reynolds number) COARE3.6-specific! 
    301          z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #LOLO: some use 1.15 not 1.1 !!! 
     301         z0t   = MIN( 1.1E-4_wp , 5.5E-5_wp*ztmp1 ) ! Scalar roughness for both theta and q (eq.28) #LB: some use 1.15 not 1.1 !!! 
    302302         z0t   = MIN( MAX(ABS(z0t), 1.E-9) , 1._wp )                      ! (prevents FPE from stupid values from masked region later on) 
    303303 
     
    315315            t_zu = t_zt - t_star/vkarmn*ztmp1 
    316316            q_zu = q_zt - q_star/vkarmn*ztmp1 
    317          END IF 
     317         ENDIF 
    318318 
    319319 
     
    330330            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    331331 
    332          END IF 
     332         ENDIF 
    333333 
    334334         IF( l_use_wl ) THEN 
     
    343343            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 
    344344            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    345          END IF 
     345         ENDIF 
    346346 
    347347         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 
    348348            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    349349            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    350          END IF 
     350         ENDIF 
    351351 
    352352      END DO !DO j_itt = 1, nb_itt 
     
    365365      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 
    366366 
    367       IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
    368       IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 
    369       IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 
    370  
    371       IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
     367      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
     368      IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 
     369      IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 
     370 
     371      IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
    372372 
    373373   END SUBROUTINE turb_coare3p0 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_coare3p6.F90

    r11962 r11963  
    6868      INTEGER :: ierr 
    6969      !!--------------------------------------------------------------------- 
    70       IF ( l_use_wl ) THEN 
     70      IF( l_use_wl ) THEN 
    7171         ierr = 0 
    7272         ALLOCATE ( Tau_ac(jpi,jpj) , Qnt_ac(jpi,jpj), dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) 
     
    7676         dT_wl(:,:)  = 0._wp 
    7777         Hz_wl(:,:)  = Hwl_max 
    78       END IF 
    79       IF ( l_use_cs ) THEN 
     78      ENDIF 
     79      IF( l_use_cs ) THEN 
    8080         ierr = 0 
    8181         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 
    8282         IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_COARE3P6_INIT => allocation of dT_cs failed!' ) 
    8383         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction 
    84       END IF 
     84      ENDIF 
    8585   END SUBROUTINE sbcblk_algo_coare3p6_init 
    8686 
     
    190190      CHARACTER(len=40), PARAMETER :: crtnm = 'turb_coare3p6@sbcblk_algo_coare3p6' 
    191191      !!---------------------------------------------------------------------------------- 
    192       IF ( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl) 
     192      IF( kt == nit000 ) CALL SBCBLK_ALGO_COARE3P6_INIT(l_use_cs, l_use_wl) 
    193193 
    194194      l_zt_equal_zu = .FALSE. 
     
    197197 
    198198      !! Initializations for cool skin and warm layer: 
    199       IF ( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
     199      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
    200200         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' ) 
    201201 
    202       IF ( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
     202      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
    203203         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' ) 
    204204 
    205       IF ( l_use_cs .OR. l_use_wl ) THEN 
     205      IF( l_use_cs .OR. l_use_wl ) THEN 
    206206         ALLOCATE ( zsst(jpi,jpj) ) 
    207207         zsst = T_s ! backing up the bulk SST 
    208208         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction 
    209209         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
    210       END IF 
     210      ENDIF 
    211211 
    212212 
     
    265265         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    266266         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    267       END IF 
     267      ENDIF 
    268268 
    269269      !! ITERATION BLOCK 
     
    288288            zeta_t = zt*ztmp0 
    289289            zeta_t = SIGN( MIN(ABS(zeta_t),50.0_wp), zeta_t ) 
    290          END IF 
     290         ENDIF 
    291291 
    292292         !! Adjustment the wind at 10m (not needed in the current algo form): 
    293          !IF ( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u)) 
     293         !IF( zu \= 10._wp ) U10 = U_zu + u_star/vkarmn*(LOG(10._wp/zu) - psi_m_coare(10._wp*ztmp0) + psi_m_coare(zeta_u)) 
    294294 
    295295         !! Roughness lengthes z0, z0t (z0q = z0t) : 
     
    315315            t_zu = t_zt - t_star/vkarmn*ztmp1 
    316316            q_zu = q_zt - q_star/vkarmn*ztmp1 
    317          END IF 
     317         ENDIF 
    318318 
    319319 
     
    330330            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    331331 
    332          END IF 
     332         ENDIF 
    333333 
    334334         IF( l_use_wl ) THEN 
     
    343343            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 
    344344            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    345          END IF 
     345         ENDIF 
    346346 
    347347         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 
    348348            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    349349            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    350          END IF 
     350         ENDIF 
    351351 
    352352      END DO !DO j_itt = 1, nb_itt 
     
    365365      IF( .NOT. l_zt_equal_zu ) DEALLOCATE( zeta_t ) 
    366366 
    367       IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
    368       IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 
    369       IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 
    370  
    371       IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
     367      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
     368      IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 
     369      IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 
     370 
     371      IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
    372372 
    373373   END SUBROUTINE turb_coare3p6 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ecmwf.F90

    r11962 r11963  
    7474      INTEGER :: ierr 
    7575      !!--------------------------------------------------------------------- 
    76       IF ( l_use_wl ) THEN 
     76      IF( l_use_wl ) THEN 
    7777         ierr = 0 
    7878         ALLOCATE ( dT_wl(jpi,jpj), Hz_wl(jpi,jpj), STAT=ierr ) 
     
    8080         dT_wl(:,:)  = 0._wp 
    8181         Hz_wl(:,:)  = rd0 ! (rd0, constant, = 3m is default for Zeng & Beljaars) 
    82       END IF 
    83       IF ( l_use_cs ) THEN 
     82      ENDIF 
     83      IF( l_use_cs ) THEN 
    8484         ierr = 0 
    8585         ALLOCATE ( dT_cs(jpi,jpj), STAT=ierr ) 
    8686         IF( ierr > 0 ) CALL ctl_stop( ' SBCBLK_ALGO_ECMWF_INIT => allocation of dT_cs failed!' ) 
    8787         dT_cs(:,:) = -0.25_wp  ! First guess of skin correction 
    88       END IF 
     88      ENDIF 
    8989   END SUBROUTINE sbcblk_algo_ecmwf_init 
    9090 
     
    195195      !!---------------------------------------------------------------------------------- 
    196196 
    197       IF ( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 
     197      IF( kt == nit000 ) CALL SBCBLK_ALGO_ECMWF_INIT(l_use_cs, l_use_wl) 
    198198 
    199199      l_zt_equal_zu = .FALSE. 
     
    201201 
    202202      !! Initializations for cool skin and warm layer: 
    203       IF ( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
     203      IF( l_use_cs .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
    204204         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use cool-skin param!' ) 
    205205 
    206       IF ( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
     206      IF( l_use_wl .AND. (.NOT.(PRESENT(Qsw) .AND. PRESENT(rad_lw) .AND. PRESENT(slp))) ) & 
    207207         &   CALL ctl_stop( '['//TRIM(crtnm)//'] => ' , 'you need to provide Qsw, rad_lw & slp to use warm-layer param!' ) 
    208208 
    209       IF ( l_use_cs .OR. l_use_wl ) THEN 
     209      IF( l_use_cs .OR. l_use_wl ) THEN 
    210210         ALLOCATE ( zsst(jpi,jpj) ) 
    211211         zsst = T_s ! backing up the bulk SST 
    212212         IF( l_use_cs ) T_s = T_s - 0.25_wp   ! First guess of correction 
    213213         q_s    = rdct_qsat_salt*q_sat(MAX(T_s, 200._wp), slp) ! First guess of q_s 
    214       END IF 
     214      ENDIF 
    215215 
    216216 
     
    270270         dt_zu = t_zu - T_s  ; dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    271271         dq_zu = q_zu - q_s  ; dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    272       END IF 
     272      ENDIF 
    273273 
    274274 
     
    293293         Linv = ztmp0*func_m*func_m/func_h / zu     ! From Eq. 3.23, Chap.3.2.3, IFS doc - Cy40r1 
    294294         !! Note: it is slightly different that the L we would get with the usual 
    295          Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) !#LOLO 
     295         Linv = SIGN( MIN(ABS(Linv),200._wp), Linv ) ! (prevent FPE from stupid values from masked region later on...) 
    296296 
    297297         !! Update func_m with new Linv: 
     
    335335            ztmp1  = LOG(zt/zu) + ztmp2 
    336336            q_zu   = q_zt - q_star/vkarmn*ztmp1 
    337          END IF 
     337         ENDIF 
    338338 
    339339         !! Updating because of updated z0 and z0t and new Linv... 
     
    355355            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    356356 
    357          END IF 
     357         ENDIF 
    358358 
    359359         IF( l_use_wl ) THEN 
     
    366366            IF( l_use_cs ) T_s(:,:) = T_s(:,:) + dT_cs(:,:)*tmask(:,:,1) 
    367367            q_s(:,:) = rdct_qsat_salt*q_sat(MAX(T_s(:,:), 200._wp), slp(:,:)) 
    368          END IF 
     368         ENDIF 
    369369 
    370370         IF( l_use_cs .OR. l_use_wl .OR. (.NOT. l_zt_equal_zu) ) THEN 
    371371            dt_zu = t_zu - T_s ;  dt_zu = SIGN( MAX(ABS(dt_zu),1.E-6_wp), dt_zu ) 
    372372            dq_zu = q_zu - q_s ;  dq_zu = SIGN( MAX(ABS(dq_zu),1.E-9_wp), dq_zu ) 
    373          END IF 
     373         ENDIF 
    374374 
    375375      END DO !DO j_itt = 1, nb_itt 
     
    384384      Cen = vkarmn*vkarmn / (log(zu/z0q)*log(zu/z0q)) 
    385385 
    386       IF ( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
    387       IF ( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 
    388       IF ( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 
    389  
    390       IF ( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
     386      IF( l_use_cs .AND. PRESENT(pdT_cs) ) pdT_cs = dT_cs 
     387      IF( l_use_wl .AND. PRESENT(pdT_wl) ) pdT_wl = dT_wl 
     388      IF( l_use_wl .AND. PRESENT(pHz_wl) ) pHz_wl = Hz_wl 
     389 
     390      IF( l_use_cs .OR. l_use_wl ) DEALLOCATE ( zsst ) 
    391391 
    392392   END SUBROUTINE turb_ecmwf 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_algo_ncar.F90

    r11845 r11963  
    171171            q_zu = q_zt - ztmp2/vkarmn*stab    ! ztmp2 is still q*      L&Y 2004 eq.(9c) 
    172172            q_zu = max(0._wp, q_zu) 
    173          END IF 
     173         ENDIF 
    174174 
    175175         ! Update neutral wind speed at 10m and neutral Cd at 10m (L&Y 2004 eq. 9a)... 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_phy.F90

    r11962 r11963  
    541541            pQns(ji,jj) = zQlat + zQsen + zQlw 
    542542 
    543             IF ( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 
     543            IF( PRESENT(Qlat) ) Qlat(ji,jj) = zQlat 
    544544         END DO 
    545545      END DO 
     
    594594            pQlat(ji,jj) = L_vap(pTs(ji,jj)) * zevap 
    595595 
    596             IF ( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 
    597             IF ( PRESENT(prhoa) ) prhoa(ji,jj) = zrho 
     596            IF( PRESENT(pEvap) ) pEvap(ji,jj) = - zevap 
     597            IF( PRESENT(prhoa) ) prhoa(ji,jj) = zrho 
    598598 
    599599         END DO 
     
    647647      pQlat = L_vap(pTs) * zevap 
    648648 
    649       IF ( PRESENT(pEvap) ) pEvap = - zevap 
    650       IF ( PRESENT(prhoa) ) prhoa = zrho 
     649      IF( PRESENT(pEvap) ) pEvap = - zevap 
     650      IF( PRESENT(prhoa) ) prhoa = zrho 
    651651 
    652652   END SUBROUTINE BULK_FORMULA_SCLR 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_coare.F90

    r11962 r11963  
    22   !!====================================================================== 
    33   !!                   ***  MODULE  sbcblk_skin_coare  *** 
    4    !! Computes: 
    5    !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer 
    6    !!     scheme used at ECMWF 
    7    !!    Using formulation/param. of COARE 3.6 (Fairall et al., 2019) 
    84   !! 
    9    !!   From routine turb_ecmwf maintained and developed in AeroBulk 
    10    !!            (https://github.com/brodeau/aerobulk) 
     5   !!   Module that gathers the cool-skin and warm-layer parameterization used 
     6   !!   in the COARE family of bulk parameterizations. 
    117   !! 
    12    !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 
     8   !!   Based on the last update for version COARE 3.6 (Fairall et al., 2019) 
     9   !! 
     10   !!   Module 'sbcblk_skin_coare' also maintained and developed in AeroBulk (as 
     11   !!            'mod_skin_coare') 
     12   !!    (https://github.com/brodeau/aerobulk)   !! 
     13   !! ** Author: L. Brodeau, November 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 
    1314   !!---------------------------------------------------------------------- 
    14    !! History :  4.0  !  2016-02  (L.Brodeau)   Original code 
    15    !!---------------------------------------------------------------------- 
    16  
    17    !!---------------------------------------------------------------------- 
    18    !!  cswl_ecmwf : computes the surface skin temperature (aka SSST) 
    19    !!               based on the cool-skin/warm-layer scheme used at ECMWF 
     15   !! History :  4.x  !  2019-11  (L.Brodeau)   Original code 
    2016   !!---------------------------------------------------------------------- 
    2117   USE oce             ! ocean dynamics and tracers 
     
    2420   USE sbc_oce         ! Surface boundary condition: ocean fields 
    2521 
    26    USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc... 
    27  
    28    USE sbcdcy          !LOLO: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE) 
    29     
     22   USE sbcblk_phy      ! misc. functions for marine ABL physics (rho_air, q_sat, bulk_formula, etc) 
     23 
     24   USE sbcdcy          !#LB: to know hour of dawn and dusk: rdawn_dcy and rdusk_dcy (needed in WL_COARE) 
     25 
    3026   USE lib_mpp         ! distribued memory computing library 
    3127   USE in_out_manager  ! I/O manager 
     
    3834 
    3935   !! Cool-skin related parameters: 
    40    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 
    41       &                        dT_cs         !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta) 
     36   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_cs !: dT due to cool-skin effect 
     37   !                                                      ! => temperature difference between air-sea interface (z=0) 
     38   !                                                      !    and right below viscous layer (z=delta) 
    4239 
    4340   !! Warm-layer related parameters: 
    44    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 
    45       &                        dT_wl,     &  !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 
    46       &                        Hz_wl,     &  !: depth of warm-layer [m] 
    47       &                        Qnt_ac,    &  !: time integral / accumulated heat stored by the warm layer Qxdt => [J/m^2] (reset to zero every midnight) 
    48       &                        Tau_ac        !: time integral / accumulated momentum Tauxdt => [N.s/m^2] (reset to zero every midnight) 
     41   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_wl !: dT due to warm-layer effect 
     42   !                                                      ! => difference between "almost surface (right below 
     43   !                                                      !    viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 
     44   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Hz_wl !: depth (aka thickness) of warm-layer [m] 
     45   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Qnt_ac !: time integral / accumulated heat stored by the warm layer 
     46   !                                                      !         Qxdt => [J/m^2] (reset to zero every midnight) 
     47   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Tau_ac !: time integral / accumulated momentum 
     48   !                                                      !         Tauxdt => [N.s/m^2] (reset to zero every midnight) 
    4949 
    5050   REAL(wp), PARAMETER, PUBLIC :: Hwl_max = 20._wp    !: maximum depth of warm layer (adjustable) 
     
    8585      !!--------------------------------------------------------------------- 
    8686      INTEGER  :: ji, jj, jc 
    87       REAL(wp) :: zQabs, zdelta, zfr 
     87      REAL(wp) :: zQabs, zdlt, zfr, zalfa, zqlat, zus 
    8888      !!--------------------------------------------------------------------- 
    8989      DO jj = 1, jpj 
     
    9191 
    9292            zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta, 
    93             !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta... 
    94  
    95             zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 
    96  
    97             DO jc = 1, 4 ! because implicit in terms of zdelta... 
    98                zfr = MAX( 0.137_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) /  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 ! 
     93            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdlt... 
     94 
     95            zalfa = alpha_sw(pSST(ji,jj)) ! (crude) thermal expansion coefficient of sea-water [1/K] 
     96            zqlat = pQlat(ji,jj) 
     97            zus   = pustar(ji,jj) 
     98 
     99 
     100            zdlt = delta_skin_layer( zalfa, zQabs, zqlat, zus ) 
     101 
     102            DO jc = 1, 4 ! because implicit in terms of zdlt... 
     103               zfr = MAX( 0.137_wp + 11._wp*zdlt & 
     104                  &       - 6.6E-5_wp/zdlt*(1._wp - EXP(-zdlt/8.E-4_wp)) & 
     105                  &      , 0.01_wp ) ! Solar absorption, Eq.16 (Fairall al. 1996b) 
     106               !                  !LB: why 0.065 and not 0.137 like in the paper??? Beljaars & Zeng use 0.065, not 0.137 ! 
    99107               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 
    100                zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pQlat(ji,jj), pustar(ji,jj) ) 
     108               zdlt = delta_skin_layer( zalfa, zQabs, zqlat, zus ) 
    101109            END DO 
    102110 
    103             dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 
     111            dT_cs(ji,jj) = zQabs*zdlt/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 
    104112 
    105113         END DO 
     
    132140      REAL(wp) :: zdTwl, zHwl, zQabs, zfr 
    133141      REAL(wp) :: zqac, ztac 
    134       REAL(wp) :: zalpha, zcd1, zcd2, flg 
     142      REAL(wp) :: zalfa, zcd1, zcd2, flg 
    135143      !!--------------------------------------------------------------------- 
    136144 
     
    142150      !! INITIALIZATION: 
    143151      zQabs  = 0._wp       ! total heat flux absorped in warm layer 
    144       zfr    = zfr0        ! initial value of solar flux absorption !LOLO: save it and use previous value !!! 
     152      zfr    = zfr0        ! initial value of solar flux absorption !#LB: save it and use previous value !!! 
    145153 
    146154      IF( .NOT. ln_dm2dc ) CALL sbc_dcy_param() ! we need to call sbc_dcy_param (sbcdcy.F90) because rdawn_dcy and rdusk_dcy are unkonwn otherwize! 
     
    148156      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 ... 
    149157 
    150       IF (lwp) THEN 
    151          WRITE(numout,*) '' 
    152          WRITE(numout,*) 'LOLO: sbcblk_skin_coare => nsec_day, ztime =', nsec_day, ztime 
    153       END IF 
    154        
    155158      DO jj = 1, jpj 
    156159         DO ji = 1, jpi 
     
    166169 
    167170            !*****  variables for warm layer  *** 
    168             zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 
    169  
    170             zcd1 = SQRT(2._wp*rich*rCp0_w/(zalpha*grav*rho0_w))        !mess-o-constants 1 
    171             zcd2 = SQRT(2._wp*zalpha*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 
    172  
    173              
     171            zalfa = alpha_sw( pSST(ji,jj) ) ! (crude) thermal expansion coefficient of sea-water [1/K] (SST accurate enough!) 
     172 
     173            zcd1 = SQRT(2._wp*rich*rCp0_w/(zalfa*grav*rho0_w))        !mess-o-constants 1 
     174            zcd2 = SQRT(2._wp*zalfa*grav/(rich*rho0_w))/(rCp0_w**1.5) !mess-o-constants 2 
     175 
     176 
    174177            znoon = MOD( 0.5_wp*(rdawn_dcy(ji,jj)+rdusk_dcy(ji,jj)), 1._wp )   ! 0<rnoon<1. => rnoon*24 = UTC time of local noon 
    175178            zmidn = MOD(              znoon-0.5_wp                 , 1._wp ) 
    176             zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight  
    177  
    178             IF ( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN 
     179            zmidn = MOD( zmidn + 0.125_wp , 1._wp ) ! 3 hours past the local midnight 
     180 
     181            IF( (ztime >= zmidn) .AND. (ztime < rdawn_dcy(ji,jj)) ) THEN 
    179182               ! Dawn reset to 0! 
    180183               l_exit       = .TRUE. 
    181184               l_destroy_wl = .TRUE. 
    182             END IF 
    183  
    184             IF ( .NOT. l_exit ) THEN 
     185            ENDIF 
     186 
     187            IF( .NOT. l_exit ) THEN 
    185188               !! Initial test on initial guess of absorbed heat flux in warm-layer: 
    186                zfr = 1._wp - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & 
    187                   &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 
    188                zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer !LOLO: depends of zfr, which is wild guess... Wrong!!! 
    189  
    190                IF ( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN 
     189               zQabs = frac_solar_abs(zHwl)*pQsw(ji,jj) + pQnsol(ji,jj) ! first guess of tot. heat flux absorbed in warm layer 
     190               !                                                        ! => #LB: depends of zfr, which is wild guess... Wrong!!! 
     191               IF( (ABS(zdTwl) < 1.E-6_wp) .AND. (zQabs <= 0._wp) ) THEN 
    191192                  ! We have not started to build a WL yet (dT==0) and there's no way it can occur now 
    192193                  ! since zQabs <= 0._wp 
    193194                  ! => no need to go further 
    194195                  l_exit = .TRUE. 
    195                END IF 
    196  
    197             END IF 
     196               ENDIF 
     197 
     198            ENDIF 
    198199 
    199200            ! Okay test on updated absorbed flux: 
    200             !LOLO: remove??? has a strong influence !!! 
    201             IF ( (.NOT.(l_exit)) .AND. (Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN 
     201            !#LB: remove??? has a strong influence !!! 
     202            IF( (.NOT. l_exit).AND.(Qnt_ac(ji,jj) + zQabs*rdt <= 0._wp) ) THEN 
    202203               l_exit       = .TRUE. 
    203204               l_destroy_wl = .TRUE. 
    204             END IF 
    205  
    206  
    207             IF ( .NOT. l_exit) THEN 
     205            ENDIF 
     206 
     207 
     208            IF( .NOT. l_exit) THEN 
    208209 
    209210               ! Two possibilities at this point: 
     
    217218               !! some part is useless if Qsw=0 !!! 
    218219               DO jl = 1, 5 
    219                   zfr = 1. - ( 0.28*0.014*(1. - EXP(-zHwl/0.014)) + 0.27*0.357*(1. - EXP(-zHwl/0.357)) & 
    220                      &        + 0.45*12.82*(1-EXP(-zHwl/12.82)) ) / zHwl 
    221                   zQabs = zfr*pQsw(ji,jj) + pQnsol(ji,jj) 
     220                  zQabs = frac_solar_abs(zHwl)*pQsw(ji,jj) + pQnsol(ji,jj) 
    222221                  zqac  = Qnt_ac(ji,jj) + zQabs*rdt ! updated heat absorbed 
    223                   IF ( zqac <= 0._wp ) EXIT 
     222                  IF( zqac <= 0._wp ) EXIT 
    224223                  zHwl = MAX( MIN( Hwl_max , zcd1*ztac/SQRT(zqac)) , 0.1_wp ) ! Warm-layer depth 
    225224               END DO 
    226225 
    227                IF ( zqac <= 0._wp ) THEN 
     226               IF( zqac <= 0._wp ) THEN 
    228227                  l_destroy_wl = .TRUE. 
    229228                  l_exit       = .TRUE. 
     
    234233                  flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl )               ! => 1 when gdept_1d(1)>zHwl (zdTwl = zdTwl) | 0 when gdept_1d(1)<zHwl (zdTwl = zdTwl*gdept_1d(1)/zHwl) 
    235234                  zdTwl = zdTwl * ( flg + (1._wp-flg)*gdept_1d(1)/zHwl ) 
    236                END IF 
    237  
    238             END IF !IF ( .NOT. l_exit) 
    239  
    240             IF ( l_destroy_wl ) THEN 
     235               ENDIF 
     236 
     237            ENDIF !IF( .NOT. l_exit) 
     238 
     239            IF( l_destroy_wl ) THEN 
    241240               zdTwl = 0._wp 
    242241               zfr   = 0.75_wp 
     
    244243               zqac  = 0._wp 
    245244               ztac  = 0._wp 
    246             END IF 
    247  
    248             IF ( iwait == 0 ) THEN 
     245            ENDIF 
     246 
     247            IF( iwait == 0 ) THEN 
    249248               !! Iteration loop within bulk algo is over, time to update what needs to be updated: 
    250249               dT_wl(ji,jj)  = zdTwl 
     
    252251               Qnt_ac(ji,jj) = zqac ! Updating Qnt_ac, heat integral 
    253252               Tau_ac(ji,jj) = ztac 
    254             END IF 
     253            ENDIF 
    255254 
    256255         END DO 
     
    276275      REAL(wp)                :: delta_skin_layer 
    277276      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!) 
    278       REAL(wp), INTENT(in)    :: pQd    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 
     277      REAL(wp), INTENT(in)    :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] 
     278      !                              !  => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 
    279279      REAL(wp), INTENT(in)    :: pQlat    ! latent heat flux [W/m^2] 
    280280      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s] 
     
    282282      REAL(wp) :: zusw, zusw2, zlamb, zQd, ztf, ztmp 
    283283      !!--------------------------------------------------------------------- 
    284        
    285       zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! LOLO: Double check sign + division by palpha !!! units are okay! 
     284 
     285      zQd = pQd + 0.026*MIN(pQlat,0._wp)*rCp0_w/rLevap/palpha   ! #LB: Double check sign + division by palpha !!! units are okay! 
    286286 
    287287      ztf = 0.5_wp + SIGN(0.5_wp, zQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case) 
     
    299299   END FUNCTION delta_skin_layer 
    300300 
     301 
     302   FUNCTION frac_solar_abs( pHwl ) 
     303      !!--------------------------------------------------------------------- 
     304      !! Fraction of solar heat flux absorbed inside warm layer 
     305      !!--------------------------------------------------------------------- 
     306      REAL(wp)             :: frac_solar_abs 
     307      REAL(wp), INTENT(in) :: pHwl   ! thickness of warm-layer [m] 
     308      !!--------------------------------------------------------------------- 
     309      frac_solar_abs = 1._wp - ( 0.28*0.014  *(1._wp - EXP(-pHwl/0.014)) & 
     310         &                       + 0.27*0.357*(1._wp - EXP(-pHwl/0.357)) & 
     311         &                       + 0.45*12.82*(1-EXP(-pHwl/12.82)) ) / pHwl 
     312   END FUNCTION frac_solar_abs 
     313 
    301314   !!====================================================================== 
    302315END MODULE sbcblk_skin_coare 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcblk_skin_ecmwf.F90

    r11826 r11963  
    22   !!====================================================================== 
    33   !!                   ***  MODULE  sbcblk_skin_ecmwf  *** 
    4    !! Computes: 
    5    !!   * the surface skin temperature (aka SSST) based on the cool-skin/warm-layer 
    6    !!     scheme used at ECMWF 
    7    !!    Using formulation/param. of IFS of ECMWF (cycle 40r1) 
    8    !!   based on IFS doc (avaible online on the ECMWF's website) 
    9    !! 
    10    !!   From routine turb_ecmwf maintained and developed in AeroBulk 
    11    !!            (https://github.com/brodeau/aerobulk) 
    12    !! 
    13    !! ** Author: L. Brodeau, September 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 
     4   !! 
     5   !!   Module that gathers the cool-skin and warm-layer parameterization used 
     6   !!   by the IFS at ECMWF (recoded from scratch => 
     7   !!   https://github.com/brodeau/aerobulk) 
     8   !! 
     9   !!  Mainly based on Zeng & Beljaars, 2005 with the more recent add-up from 
     10   !!  Takaya et al., 2010 when it comes to the warm-layer parameterization 
     11   !!  (contribution of extra mixing due to Langmuir circulation) 
     12   !! 
     13   !!  - Zeng X., and A. Beljaars, 2005: A prognostic scheme of sea surface skin 
     14   !!    temperature for modeling and data assimilation. Geophysical Research 
     15   !!    Letters, 32 (14) , pp. 1-4. 
     16   !! 
     17   !!  - Takaya, Y., J.-R. Bildot, A. C. M. Beljaars, and P. A. E. M. Janssen, 
     18   !!    2010: Refinements to a prognostic scheme of skin sea surface 
     19   !!    temperature. J. Geophys. Res., 115, C06009, doi:10.1029/2009JC005985 
     20   !! 
     21   !!   Most of the formula are taken from the documentation of IFS of ECMWF 
     22   !!            (cycle 40r1) (avaible online on the ECMWF's website) 
     23   !! 
     24   !!   Routine 'sbcblk_skin_ecmwf' also maintained and developed in AeroBulk (as 
     25   !!            'mod_skin_ecmwf') 
     26   !!    (https://github.com/brodeau/aerobulk) 
     27   !! 
     28   !! ** Author: L. Brodeau, November 2019 / AeroBulk (https://github.com/brodeau/aerobulk) 
    1429   !!---------------------------------------------------------------------- 
    15    !! History :  4.0  !  2016-02  (L.Brodeau)   Original code 
    16    !!---------------------------------------------------------------------- 
    17  
    18    !!---------------------------------------------------------------------- 
    19    !!  cswl_ecmwf : computes the surface skin temperature (aka SSST) 
    20    !!               based on the cool-skin/warm-layer scheme used at ECMWF 
     30   !! History :  4.x  !  2019-11  (L.Brodeau)   Original code 
    2131   !!---------------------------------------------------------------------- 
    2232   USE oce             ! ocean dynamics and tracers 
     
    2535   USE sbc_oce         ! Surface boundary condition: ocean fields 
    2636 
    27    USE sbcblk_phy      !LOLO: all thermodynamics functions, rho_air, q_sat, etc... 
     37   USE sbcblk_phy      ! misc. functions for marine ABL physics (rho_air, q_sat, bulk_formula, etc) 
    2838 
    2939   USE lib_mpp         ! distribued memory computing library 
     
    3141   USE lib_fortran     ! to use key_nosignedzero 
    3242 
    33  
    3443   IMPLICIT NONE 
    3544   PRIVATE 
     
    3847 
    3948   !! Cool-skin related parameters: 
    40    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 
    41       &                        dT_cs         !: dT due to cool-skin effect => temperature difference between air-sea interface (z=0) and right below viscous layer (z=delta) 
     49   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_cs !: dT due to cool-skin effect 
     50   !                                                      ! => temperature difference between air-sea interface (z=0) 
     51   !                                                      !    and right below viscous layer (z=delta) 
    4252 
    4353   !! Warm-layer related parameters: 
    44    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: & 
    45       &                        dT_wl,  &   !: dT due to warm-layer effect => difference between "almost surface (right below viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 
    46       &                        Hz_wl       !: depth of warm-layer [m] 
     54   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: dT_wl !: dT due to warm-layer effect 
     55   !                                                      ! => difference between "almost surface (right below 
     56   !                                                      !    viscous layer, z=delta) and depth of bulk SST (z=gdept_1d(1)) 
     57   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:), PUBLIC :: Hz_wl !: depth (aka thickness) of warm-layer [m] 
    4758   ! 
    4859   REAL(wp), PARAMETER, PUBLIC :: rd0  = 3.    !: Depth scale [m] of warm layer, "d" in Eq.11 (Zeng & Beljaars 2005) 
     
    6273      !! Cool-skin parameterization, based on Fairall et al., 1996: 
    6374      !! 
    64       !! Fairall, C. W., Bradley, E. F., Godfrey, J. S., Wick, G. A., 
    65       !! Edson, J. B., and Young, G. S. ( 1996), Cool‐skin and warm‐layer 
    66       !! effects on sea surface temperature, J. Geophys. Res., 101( C1), 1295-1308, 
    67       !! doi:10.1029/95JC03190. 
     75      !!  - Zeng X., and A. Beljaars, 2005: A prognostic scheme of sea surface 
     76      !!    skin temperature for modeling and data assimilation. Geophysical 
     77      !!    Research Letters, 32 (14) , pp. 1-4. 
    6878      !! 
    6979      !!------------------------------------------------------------------ 
     
    8191      !!--------------------------------------------------------------------- 
    8292      INTEGER  :: ji, jj, jc 
    83       REAL(wp) :: zQabs, zdelta, zfr 
     93      REAL(wp) :: zQabs, zdlt, zfr, zalfa, zus 
    8494      !!--------------------------------------------------------------------- 
    8595      DO jj = 1, jpj 
     
    8797 
    8898            zQabs = pQnsol(ji,jj) ! first guess of heat flux absorbed within the viscous sublayer of thicknes delta, 
    89             !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdelta... 
    90  
    91             zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 
    92  
    93             DO jc = 1, 4 ! because implicit in terms of zdelta... 
    94                zfr = MAX( 0.065_wp + 11._wp*zdelta - 6.6E-5_wp/zdelta*(1._wp - EXP(-zdelta/8.E-4_wp)) , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005 
    95                !              =>  (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996) 
     99            !                     !   => we DO not miss a lot assuming 0 solar flux absorbed in the tiny layer of thicknes zdlt... 
     100 
     101            zalfa = alpha_sw(pSST(ji,jj)) ! (crude) thermal expansion coefficient of sea-water [1/K] 
     102            zus   = pustar(ji,jj) 
     103 
     104            zdlt = delta_skin_layer( zalfa, zQabs, zus ) 
     105 
     106            DO jc = 1, 4 ! because implicit in terms of zdlt... 
     107               zfr = MAX( 0.065_wp + 11._wp*zdlt  & 
     108                  &       - 6.6E-5_wp/zdlt*(1._wp - EXP(-zdlt/8.E-4_wp)) & 
     109                  &      , 0.01_wp ) ! Solar absorption, Eq.(5) Zeng & Beljaars, 2005 
     110               !                     !  =>  (WARNING: 0.065 rather than 0.137 in Fairal et al. 1996) 
    96111               zQabs = pQnsol(ji,jj) + zfr*pQsw(ji,jj) 
    97                zdelta = delta_skin_layer( alpha_sw(pSST(ji,jj)), zQabs, pustar(ji,jj) ) 
     112               zdlt = delta_skin_layer( zalfa, zQabs, zus ) 
    98113            END DO 
    99114 
    100             dT_cs(ji,jj) = zQabs*zdelta/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 
     115            dT_cs(ji,jj) = zQabs*zdlt/rk0_w   ! temperature increment, yes dT_cs can actually > 0, if Qabs > 0 (rare but possible!) 
    101116 
    102117         END DO 
     
    106121 
    107122 
    108  
    109123   SUBROUTINE WL_ECMWF( pQsw, pQnsol, pustar, pSST,  pustk ) 
    110124      !!--------------------------------------------------------------------- 
    111125      !! 
    112       !!  Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) 
    113       !!  " A prognostic scheme of sea surface skin temperature for modeling and data assimilation " 
     126      !!  Warm-Layer scheme according to Zeng & Beljaars, 2005 (GRL) with the 
     127      !!  more recent add-up from Takaya et al., 2010 when it comes to the 
     128      !!  warm-layer parameterization (contribution of extra mixing due to 
     129      !!  Langmuir circulation) 
     130      !! 
     131      !!  - Zeng X., and A. Beljaars, 2005: A prognostic scheme of sea surface skin 
     132      !!    temperature for modeling and data assimilation. Geophysical Research 
     133      !!    Letters, 32 (14) , pp. 1-4. 
     134      !! 
     135      !!  - Takaya, Y., J.-R. Bildot, A. C. M. Beljaars, and P. A. E. M. Janssen, 
     136      !!    2010: Refinements to a prognostic scheme of skin sea surface 
     137      !!    temperature. J. Geophys. Res., 115, C06009, doi:10.1029/2009JC005985 
    114138      !! 
    115139      !!  STIL NO PROGNOSTIC EQUATION FOR THE DEPTH OF THE WARM-LAYER! 
    116140      !! 
    117       !!    As included in IFS Cy45r1   /  E.C.M.W.F. 
    118141      !!     ------------------------------------------------------------------ 
    119142      !! 
     
    133156      INTEGER :: ji, jj, jc 
    134157      ! 
    135       REAL(wp) :: & 
    136          & zHwl,    &  !: thickness of the warm-layer [m] 
    137          & ztcorr,  &  !: correction of dT w.r.t measurement depth of bulk SST (first T-point) 
    138          & zalpha, & !: thermal expansion coefficient of sea-water [1/K] 
    139          & zdTwl_b, zdTwl_n, & ! temp. diff. between "almost surface (right below viscous layer) and bottom of WL 
    140          & zfr, zeta, & 
    141          & zusw, zusw2, & 
    142          & zLa, zfLa, & 
    143          & flg, zwf, zQabs, & 
    144          & ZA, ZB, zL1, zL2, & 
    145          &  zcst0, zcst1, zcst2, zcst3 
     158      REAL(wp) :: zHwl      !: thickness of the warm-layer [m] 
     159      REAL(wp) :: ztcorr    !: correction of dT w.r.t measurement depth of bulk SST (first T-point) 
     160      REAL(wp) :: zalfa     !: thermal expansion coefficient of sea-water [1/K] 
     161      REAL(wp) :: zdTwl_b, zdTwl_n  !: temp. diff. between "almost surface (right below viscous layer) and bottom of WL 
     162      REAL(wp) :: zfr, zeta  
     163      REAL(wp) :: zusw, zusw2  
     164      REAL(wp) :: zLa, zfLa  
     165      REAL(wp) :: flg, zwf, zQabs  
     166      REAL(wp) :: ZA, ZB, zL1, zL2 
     167      REAL(wp) :: zcst0, zcst1, zcst2, zcst3 
    146168      ! 
    147169      LOGICAL :: l_pustk_known 
     
    149171 
    150172      l_pustk_known = .FALSE. 
    151       IF ( PRESENT(pustk) ) l_pustk_known = .TRUE. 
     173      IF( PRESENT(pustk) ) l_pustk_known = .TRUE. 
    152174 
    153175      DO jj = 1, jpj 
     
    158180 
    159181            !! Previous value of dT / warm-layer, adapted to depth: 
    160             flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl )               ! => 1 when gdept_1d(1)>zHwl (dT_wl(ji,jj) = zdTwl) | 0 when z_s$ 
     182            flg = 0.5_wp + SIGN( 0.5_wp , gdept_1d(1)-zHwl ) ! => 1 when gdept_1d(1)>zHwl (dT_wl(ji,jj) = zdTwl) | 0 when z_s$ 
    161183            ztcorr = flg + (1._wp - flg)*gdept_1d(1)/zHwl 
    162184            zdTwl_b = MAX ( dT_wl(ji,jj) / ztcorr , 0._wp ) 
     
    166188            !! => so here since pdT is difference between surface and gdept_1d(1), need to increase fof zdTwl ! 
    167189 
    168             zalpha = alpha_sw( pSST(ji,jj) ) ! thermal expansion coefficient of sea-water (SST accurate enough!) 
    169  
     190            zalfa = alpha_sw( pSST(ji,jj) ) ! (crude) thermal expansion coefficient of sea-water [1/K] (SST accurate enough!) 
    170191 
    171192            ! *** zfr = Fraction of solar radiation absorbed in warm layer (-) 
     
    178199 
    179200            ! Langmuir: 
    180             IF ( l_pustk_known ) THEN 
     201            IF( l_pustk_known ) THEN 
    181202               zLa = SQRT(zusw/MAX(pustk(ji,jj),1.E-6)) 
    182203            ELSE 
    183204               zla = 0.3_wp 
    184             END IF 
     205            ENDIF 
    185206            zfLa = MAX( zla**(-2._wp/3._wp) , 1._wp )   ! Eq.(6) 
    186207 
    187208            zwf = 0.5_wp + SIGN(0.5_wp, zQabs)  ! zQabs > 0. => 1.  / zQabs < 0. => 0. 
    188209 
    189             zcst1 = vkarmn*grav*zalpha 
     210            zcst1 = vkarmn*grav*zalfa 
    190211 
    191212            ! 1/L when zQabs > 0 : 
    192213            zL2 = zcst1*zQabs / (zRhoCp_w*zusw2*zusw) 
    193                 
     214 
    194215            zcst2 = zcst1 / ( 5._wp*zHwl*zusw2 )  !OR: zcst2 = zcst1*rNuwl0 / ( 5._wp*zHwl*zusw2 ) ??? 
    195216 
    196217            zcst0 = rdt * (rNuwl0 + 1._wp) / zHwl 
    197              
     218 
    198219            ZA = zcst0 * zQabs / ( rNuwl0 * zRhoCp_w ) 
    199220 
     
    203224            !! the sake of optimizations... So all these operations are not done 
    204225            !! over and over within the iteration loop... 
    205              
     226 
    206227            !! T R U L L Y   I M P L I C I T => needs itteration 
    207228            !! => have to itterate just because the 1/(Monin-Obukhov length), zL1, uses zdTwl when zQabs < 0.. 
     
    209230            zdTwl_n = zdTwl_b 
    210231            DO jc = 1, 10 
    211                 
     232 
    212233               zdTwl_n = 0.5_wp * ( zdTwl_n + zdTwl_b ) ! semi implicit, for faster convergence 
    213                 
     234 
    214235               ! 1/L when zdTwl > 0 .AND. zQabs < 0 : 
    215                zL1 =         SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalpha/( 5._wp*zHwl ) ) / zusw 
    216                !zL1 = vkarmn*SQRT( zdTwl_n       *grav*zalpha        / ( 5._wp*zHwl ) ) / zusw   ! => vkarmn outside, not inside zcst1 (just for this particular line) ??? 
    217                 
     236               zL1 =         SQRT( zdTwl_n * zcst2 ) ! / zusw !!! Or??? => vkarmn * SQRT( zdTwl_n*grav*zalfa/( 5._wp*zHwl ) ) / zusw 
     237 
    218238               ! Stability parameter (z/L): 
    219239               zeta =  (1._wp - zwf) * zHwl*zL1   +   zwf * zHwl*zL2 
     
    224244 
    225245            END DO 
    226              
     246 
    227247            !! Update: 
    228248            dT_wl(ji,jj) = zdTwl_n * ztcorr 
    229              
     249 
    230250         END DO 
    231251      END DO 
    232252 
    233253   END SUBROUTINE WL_ECMWF 
    234  
    235254 
    236255 
     
    249268      REAL(wp)                :: delta_skin_layer 
    250269      REAL(wp), INTENT(in)    :: palpha   ! thermal expansion coefficient of sea-water (SST accurate enough!) 
    251       REAL(wp), INTENT(in)    :: pQd    ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 
     270      REAL(wp), INTENT(in)    :: pQd ! < 0 !!! part of the net heat flux actually absorbed in the WL [W/m^2] 
     271      !                              !  => term "Q + Rs*fs" in eq.6 of Fairall et al. 1996 
    252272      REAL(wp), INTENT(in)    :: pustar_a ! friction velocity in the air (u*) [m/s] 
    253273      !!--------------------------------------------------------------------- 
     
    255275      !!--------------------------------------------------------------------- 
    256276      ztf = 0.5_wp + SIGN(0.5_wp, pQd)  ! Qabs < 0 => cooling of the viscous layer => ztf = 0 (regular case) 
    257       !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 (ex: weak evaporation and strong positive sensible heat flux) 
    258       ! 
     277      !                                 ! Qabs > 0 => warming of the viscous layer => ztf = 1 
     278      !                                 !    (ex: weak evaporation and strong positive sensible heat flux) 
    259279      zusw  = MAX(pustar_a, 1.E-4_wp) * sq_radrw    ! u* in the water 
    260280      zusw2 = zusw*zusw 
     
    281301      REAL(wp) :: ztf, zzt2 
    282302      !!--------------------------------------------------------------------- 
    283       ! 
    284303      zzt2 = pzeta*pzeta 
    285       ! 
    286304      ztf = 0.5_wp + SIGN(0.5_wp, pzeta)  ! zeta > 0 => ztf = 1 
    287305      !                                   ! zeta < 0 => ztf = 0 
    288306      PHI =      ztf     * ( 1. + (5.*pzeta + 4.*zzt2)/(1. + 3.*pzeta + 0.25*zzt2) ) &   ! zeta > 0 
    289307         &  + (1. - ztf) * 1./SQRT( 1. - 16.*(-ABS(pzeta)) )                             ! zeta < 0 
    290       ! 
    291308   END FUNCTION PHI 
    292309 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbccpl.F90

    r11831 r11963  
    453453      CASE( 'conservative'  ) 
    454454         srcv( (/jpr_rain, jpr_snow, jpr_ievp, jpr_tevp/) )%laction = .TRUE. 
    455          IF ( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE. 
     455         IF( k_ice <= 1 )  srcv(jpr_ievp)%laction = .FALSE. 
    456456      CASE( 'oce and ice'   )   ;   srcv( (/jpr_ievp, jpr_sbpr, jpr_semp, jpr_oemp/) )%laction = .TRUE. 
    457457      CASE default              ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_rcv_emp%cldes' ) 
     
    558558      srcv(jpr_botm )%clname = 'OBotMlt' 
    559559      IF( TRIM(sn_rcv_iceflx%cldes) == 'coupled' ) THEN 
    560          IF ( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN 
     560         IF( TRIM( sn_rcv_iceflx%clcat ) == 'yes' ) THEN 
    561561            srcv(jpr_topm:jpr_botm)%nct = nn_cats_cpl 
    562562         ELSE 
     
    569569      !                                                      ! ------------------------- ! 
    570570      srcv(jpr_ts_ice)%clname = 'OTsfIce'    ! needed by Met Office 
    571       IF ( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' )   srcv(jpr_ts_ice)%laction = .TRUE. 
    572       IF ( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' )   srcv(jpr_ts_ice)%nct     = nn_cats_cpl 
    573       IF ( TRIM( sn_rcv_emp%clcat    ) == 'yes' )   srcv(jpr_ievp)%nct       = nn_cats_cpl 
     571      IF( TRIM( sn_rcv_ts_ice%cldes ) == 'ice' )   srcv(jpr_ts_ice)%laction = .TRUE. 
     572      IF( TRIM( sn_rcv_ts_ice%clcat ) == 'yes' )   srcv(jpr_ts_ice)%nct     = nn_cats_cpl 
     573      IF( TRIM( sn_rcv_emp%clcat    ) == 'yes' )   srcv(jpr_ievp)%nct       = nn_cats_cpl 
    574574 
    575575      !                                                      ! ------------------------- ! 
     
    693693         ! for example O_Runoff received by OPA from SAS and therefore O_Runoff received by SAS from the Atmosphere 
    694694         DO jn = 1, jprcv 
    695             IF ( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname)) 
     695            IF( srcv(jn)%clname(1:1) == "O" ) srcv(jn)%clname = "S"//srcv(jn)%clname(2:LEN(srcv(jn)%clname)) 
    696696         END DO 
    697697         ! 
     
    720720      ! =================================================== ! 
    721721      DO jn = 1, jprcv 
    722          IF ( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) ) 
     722         IF( srcv(jn)%laction ) ALLOCATE( frcv(jn)%z3(jpi,jpj,srcv(jn)%nct) ) 
    723723      END DO 
    724724      ! Allocate taum part of frcv which is used even when not received as coupling field 
    725       IF ( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) ) 
     725      IF( .NOT. srcv(jpr_taum)%laction ) ALLOCATE( frcv(jpr_taum)%z3(jpi,jpj,srcv(jpr_taum)%nct) ) 
    726726      ! Allocate w10m part of frcv which is used even when not received as coupling field 
    727       IF ( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) ) 
     727      IF( .NOT. srcv(jpr_w10m)%laction ) ALLOCATE( frcv(jpr_w10m)%z3(jpi,jpj,srcv(jpr_w10m)%nct) ) 
    728728      ! Allocate jpr_otx1 part of frcv which is used even when not received as coupling field 
    729       IF ( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) ) 
    730       IF ( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) ) 
     729      IF( .NOT. srcv(jpr_otx1)%laction ) ALLOCATE( frcv(jpr_otx1)%z3(jpi,jpj,srcv(jpr_otx1)%nct) ) 
     730      IF( .NOT. srcv(jpr_oty1)%laction ) ALLOCATE( frcv(jpr_oty1)%z3(jpi,jpj,srcv(jpr_oty1)%nct) ) 
    731731      ! Allocate itx1 and ity1 as they are used in sbc_cpl_ice_tau even if srcv(jpr_itx1)%laction = .FALSE. 
    732732      IF( k_ice /= 0 ) THEN 
    733          IF ( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) ) 
    734          IF ( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) ) 
    735       END IF 
     733         IF( .NOT. srcv(jpr_itx1)%laction ) ALLOCATE( frcv(jpr_itx1)%z3(jpi,jpj,srcv(jpr_itx1)%nct) ) 
     734         IF( .NOT. srcv(jpr_ity1)%laction ) ALLOCATE( frcv(jpr_ity1)%z3(jpi,jpj,srcv(jpr_ity1)%nct) ) 
     735      ENDIF 
    736736 
    737737      ! ================================ ! 
     
    757757      CASE( 'oce and ice' , 'weighted oce and ice' , 'oce and weighted ice' ) 
    758758         ssnd( (/jps_toce, jps_tice/) )%laction = .TRUE. 
    759          IF ( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = nn_cats_cpl 
     759         IF( TRIM( sn_snd_temp%clcat ) == 'yes' )  ssnd(jps_tice)%nct = nn_cats_cpl 
    760760      CASE( 'mixed oce-ice'                        )   ;   ssnd( jps_tmix )%laction = .TRUE. 
    761761      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_temp%cldes' ) 
     
    777777      !     1. sending mixed oce-ice albedo or 
    778778      !     2. receiving mixed oce-ice solar radiation  
    779       IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN 
     779      IF( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN 
    780780         CALL oce_alb( zaos, zacs ) 
    781781         ! Due to lack of information on nebulosity : mean clear/overcast sky 
     
    796796         ssnd(jps_fice1)%laction = .TRUE.                 ! First-order regridded ice concentration, to be used producing atmos-to-ice fluxes (Met Office requirement) 
    797797! Currently no namelist entry to determine sending of multi-category ice fraction so use the thickness entry for now 
    798          IF ( TRIM( sn_snd_thick%clcat  ) == 'yes' ) ssnd(jps_fice)%nct  = nn_cats_cpl 
    799          IF ( TRIM( sn_snd_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = nn_cats_cpl 
     798         IF( TRIM( sn_snd_thick%clcat  ) == 'yes' ) ssnd(jps_fice)%nct  = nn_cats_cpl 
     799         IF( TRIM( sn_snd_thick1%clcat ) == 'yes' ) ssnd(jps_fice1)%nct = nn_cats_cpl 
    800800      ENDIF 
    801801       
    802       IF (TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE.  
     802      IF(TRIM( sn_snd_ifrac%cldes )  == 'coupled') ssnd(jps_ficet)%laction = .TRUE.  
    803803 
    804804      SELECT CASE ( TRIM( sn_snd_thick%cldes ) ) 
     
    806806      CASE( 'ice and snow' )  
    807807         ssnd(jps_hice:jps_hsnw)%laction = .TRUE. 
    808          IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN 
     808         IF( TRIM( sn_snd_thick%clcat ) == 'yes' ) THEN 
    809809            ssnd(jps_hice:jps_hsnw)%nct = nn_cats_cpl 
    810810         ENDIF 
    811811      CASE ( 'weighted ice and snow' )  
    812812         ssnd(jps_hice:jps_hsnw)%laction = .TRUE. 
    813          IF ( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = nn_cats_cpl 
     813         IF( TRIM( sn_snd_thick%clcat ) == 'yes' ) ssnd(jps_hice:jps_hsnw)%nct = nn_cats_cpl 
    814814      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_thick%cldes' ) 
    815815      END SELECT 
     
    828828         ssnd(jps_a_p)%laction  = .TRUE.  
    829829         ssnd(jps_ht_p)%laction = .TRUE.  
    830          IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN  
     830         IF( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN  
    831831            ssnd(jps_a_p)%nct  = nn_cats_cpl  
    832832            ssnd(jps_ht_p)%nct = nn_cats_cpl  
    833833         ELSE  
    834             IF ( nn_cats_cpl > 1 ) THEN  
     834            IF( nn_cats_cpl > 1 ) THEN  
    835835               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_mpnd%cldes if not exchanging category fields' )  
    836836            ENDIF  
     
    839839         ssnd(jps_a_p)%laction  = .TRUE.  
    840840         ssnd(jps_ht_p)%laction = .TRUE.  
    841          IF ( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN  
     841         IF( TRIM( sn_snd_mpnd%clcat ) == 'yes' ) THEN  
    842842            ssnd(jps_a_p)%nct  = nn_cats_cpl   
    843843            ssnd(jps_ht_p)%nct = nn_cats_cpl   
     
    914914      CASE ( 'ice only' )  
    915915         ssnd(jps_ttilyr)%laction = .TRUE.  
    916          IF ( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) THEN  
     916         IF( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) THEN  
    917917            ssnd(jps_ttilyr)%nct = nn_cats_cpl  
    918918         ELSE  
    919             IF ( nn_cats_cpl > 1 ) THEN  
     919            IF( nn_cats_cpl > 1 ) THEN  
    920920               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_ttilyr%cldes if not exchanging category fields' )  
    921921            ENDIF  
     
    923923      CASE ( 'weighted ice' )  
    924924         ssnd(jps_ttilyr)%laction = .TRUE.  
    925          IF ( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) ssnd(jps_ttilyr)%nct = nn_cats_cpl  
     925         IF( TRIM( sn_snd_ttilyr%clcat ) == 'yes' ) ssnd(jps_ttilyr)%nct = nn_cats_cpl  
    926926      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_ttilyr%cldes;'//sn_snd_ttilyr%cldes )  
    927927      END SELECT  
     
    933933      CASE ( 'ice only' )  
    934934         ssnd(jps_kice)%laction = .TRUE.  
    935          IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN  
     935         IF( TRIM( sn_snd_cond%clcat ) == 'yes' ) THEN  
    936936            ssnd(jps_kice)%nct = nn_cats_cpl  
    937937         ELSE  
    938             IF ( nn_cats_cpl > 1 ) THEN  
     938            IF( nn_cats_cpl > 1 ) THEN  
    939939               CALL ctl_stop( 'sbc_cpl_init: use weighted ice option for sn_snd_cond%cldes if not exchanging category fields' )  
    940940            ENDIF  
     
    942942      CASE ( 'weighted ice' )  
    943943         ssnd(jps_kice)%laction = .TRUE.  
    944          IF ( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = nn_cats_cpl  
     944         IF( TRIM( sn_snd_cond%clcat ) == 'yes' ) ssnd(jps_kice)%nct = nn_cats_cpl  
    945945      CASE default   ;   CALL ctl_stop( 'sbc_cpl_init: wrong definition of sn_snd_cond%cldes;'//sn_snd_cond%cldes )  
    946946      END SELECT  
     
    10031003         ! for example O_SSTSST sent by OPA to SAS and therefore S_SSTSST sent by SAS to the Atmosphere 
    10041004         DO jn = 1, jpsnd 
    1005             IF ( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname)) 
     1005            IF( ssnd(jn)%clname(1:1) == "O" ) ssnd(jn)%clname = "S"//ssnd(jn)%clname(2:LEN(ssnd(jn)%clname)) 
    10061006         END DO 
    10071007         ! 
     
    10301030      CALL cpl_define(jprcv, jpsnd, nn_cplmodel) 
    10311031       
    1032       IF (ln_usecplmask) THEN  
     1032      IF(ln_usecplmask) THEN  
    10331033         xcplmask(:,:,:) = 0. 
    10341034         CALL iom_open( 'cplmask', inum ) 
     
    12661266     
    12671267          IF( kt == nit000 ) ssh_ibb(:,:) = ssh_ib(:,:)  ! correct this later (read from restart if possible)  
    1268       END IF  
     1268      ENDIF  
    12691269      ! 
    12701270      IF( ln_sdw ) THEN  ! Stokes Drift correction activated 
     
    14151415         ELSE IF( srcv(jpr_qnsmix)%laction ) THEN   ;   zqns(:,:) = frcv(jpr_qnsmix)%z3(:,:,1) 
    14161416         ELSE                                       ;   zqns(:,:) = 0._wp 
    1417          END IF 
     1417         ENDIF 
    14181418         ! update qns over the free ocean with: 
    14191419         IF( nn_components /= jp_iam_opa ) THEN 
     
    16871687      ! --- evaporation over ice (kg/m2/s) --- ! 
    16881688      DO jl=1,jpl 
    1689          IF (sn_rcv_emp%clcat == 'yes') THEN   ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl) 
     1689         IF(sn_rcv_emp%clcat == 'yes') THEN   ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,jl) 
    16901690         ELSE                                  ;   zevap_ice(:,:,jl) = frcv(jpr_ievp)%z3(:,:,1 )   ;   ENDIF 
    16911691      ENDDO 
     
    17861786      CASE( 'conservative' )     ! the required fields are directly provided 
    17871787         zqns_tot(:,:) = frcv(jpr_qnsmix)%z3(:,:,1) 
    1788          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN 
     1788         IF( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN 
    17891789            zqns_ice(:,:,1:jpl) = frcv(jpr_qnsice)%z3(:,:,1:jpl) 
    17901790         ELSE 
     
    17951795      CASE( 'oce and ice' )      ! the total flux is computed from ocean and ice fluxes 
    17961796         zqns_tot(:,:) =  ziceld(:,:) * frcv(jpr_qnsoce)%z3(:,:,1) 
    1797          IF ( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN 
     1797         IF( TRIM(sn_rcv_qns%clcat) == 'yes' ) THEN 
    17981798            DO jl=1,jpl 
    17991799               zqns_tot(:,:   ) = zqns_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qnsice)%z3(:,:,jl)    
     
    18971897#endif 
    18981898      ! outputs 
    1899       IF ( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )                      ! latent heat from calving 
    1900       IF ( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting 
    1901       IF ( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average) 
    1902       IF ( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) & 
     1899      IF( srcv(jpr_cal)%laction       ) CALL iom_put('hflx_cal_cea'    , - frcv(jpr_cal)%z3(:,:,1) * rLfus )                      ! latent heat from calving 
     1900      IF( srcv(jpr_icb)%laction       ) CALL iom_put('hflx_icb_cea'    , - frcv(jpr_icb)%z3(:,:,1) * rLfus )                      ! latent heat from icebergs melting 
     1901      IF( iom_use('hflx_rain_cea')    ) CALL iom_put('hflx_rain_cea'   , ( tprecip(:,:) - sprecip(:,:) ) * zcptrain(:,:) )        ! heat flux from rain (cell average) 
     1902      IF( iom_use('hflx_evap_cea')    ) CALL iom_put('hflx_evap_cea'   , ( frcv(jpr_tevp)%z3(:,:,1) - frcv(jpr_ievp)%z3(:,:,1) & 
    19031903           &                                                              * picefr(:,:) ) * zcptn(:,:) * tmask(:,:,1) )            ! heat flux from evap (cell average) 
    1904       IF ( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average) 
    1905       IF ( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 
     1904      IF( iom_use('hflx_snow_cea')    ) CALL iom_put('hflx_snow_cea'   , sprecip(:,:) * ( zcptsnw(:,:) - rLfus )  )               ! heat flux from snow (cell average) 
     1905      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) & 
    19061906           &                                                              * ( 1._wp - zsnw(:,:) )                  )               ! heat flux from snow (over ocean) 
    1907       IF ( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &  
     1907      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea', sprecip(:,:) * ( zcptsnw(:,:) - rLfus ) &  
    19081908           &                                                              *           zsnw(:,:)                    )               ! heat flux from snow (over ice) 
    19091909      ! note: hflx for runoff and iceshelf are done in sbcrnf and sbcisf resp. 
     
    19161916      CASE( 'conservative' ) 
    19171917         zqsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1) 
    1918          IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN 
     1918         IF( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN 
    19191919            zqsr_ice(:,:,1:jpl) = frcv(jpr_qsrice)%z3(:,:,1:jpl) 
    19201920         ELSE 
     
    19281928      CASE( 'oce and ice' ) 
    19291929         zqsr_tot(:,:  ) =  ziceld(:,:) * frcv(jpr_qsroce)%z3(:,:,1) 
    1930          IF ( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN 
     1930         IF( TRIM(sn_rcv_qsr%clcat) == 'yes' ) THEN 
    19311931            DO jl = 1, jpl 
    19321932               zqsr_tot(:,:   ) = zqsr_tot(:,:) + a_i(:,:,jl) * frcv(jpr_qsrice)%z3(:,:,jl)    
     
    19841984      !                                                      ! ========================= ! 
    19851985      CASE ('coupled') 
    1986          IF ( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN 
     1986         IF( TRIM(sn_rcv_dqnsdt%clcat) == 'yes' ) THEN 
    19871987            zdqns_ice(:,:,1:jpl) = frcv(jpr_dqnsdt)%z3(:,:,1:jpl) 
    19881988         ELSE 
     
    20622062      IF( ssnd(jps_toce)%laction .OR. ssnd(jps_tice)%laction .OR. ssnd(jps_tmix)%laction ) THEN 
    20632063          
    2064          IF ( nn_components == jp_iam_opa ) THEN 
     2064         IF( nn_components == jp_iam_opa ) THEN 
    20652065            ztmp1(:,:) = tsn(:,:,1,jp_tem)   ! send temperature as it is (potential or conservative) -> use of l_useCT on the received part 
    20662066         ELSE 
     
    24672467      IF( ssnd(jps_ficet)%laction ) THEN  
    24682468         CALL cpl_snd( jps_ficet, isec, RESHAPE ( fr_i, (/jpi,jpj,1/) ), info )  
    2469       END IF  
     2469      ENDIF  
    24702470      !                                                      ! ------------------------- !  
    24712471      !                                                      !   Water levels to waves   !  
     
    24822482         ENDIF   
    24832483         CALL cpl_snd( jps_wlev  , isec, RESHAPE ( ztmp1, (/jpi,jpj,1/) ), info )  
    2484       END IF  
     2484      ENDIF  
    24852485      ! 
    24862486      !  Fields sent by OPA to SAS when doing OPA<->SAS coupling 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcdcy.F90

    r11892 r11963  
    9494            WRITE(numout,*) 
    9595         ENDIF 
    96       END IF 
     96      ENDIF 
    9797 
    9898      ! Setting parameters for each new day: 
     
    121121                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
    122122                  ztmpm = zupusd - zlousd 
    123                   IF ( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 
     123                  IF( ztmpm .EQ. 0 ) imask_night(ji,jj) = 1 
    124124                  ! 
    125125               ELSE                                         ! day time in two parts 
     
    135135                  ztmpm = ztmpm1 + ztmpm2 
    136136                  zqsrout(ji,jj) = pqsrin(ji,jj) * ztmp * rscal(ji,jj) 
    137                   IF (ztmpm .EQ. 0.) imask_night(ji,jj) = 1 
     137                  IF(ztmpm .EQ. 0.) imask_night(ji,jj) = 1 
    138138               ENDIF 
    139139            ELSE                                   ! 24h light or 24h night 
     
    206206         DO jj = 1, jpj 
    207207            DO ji = 1, jpi 
    208                IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
     208               IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
    209209                  ! When is it night? 
    210210                  ztx = 1._wp/(2._wp*rpi) * (ACOS(rab(ji,jj)) - rcc(ji,jj)) 
    211211                  ztest = -rbb(ji,jj) * SIN( rcc(ji,jj) + 2._wp*rpi * ztx ) 
    212212                  ! is it dawn or dusk? 
    213                   IF ( ztest > 0._wp ) THEN 
     213                  IF( ztest > 0._wp ) THEN 
    214214                     rdawn_dcy(ji,jj) = ztx 
    215215                     rdusk_dcy(ji,jj) = rtmd(ji,jj) + ( rtmd(ji,jj) - rdawn_dcy(ji,jj) ) 
     
    232232         DO jj = 1, jpj 
    233233            DO ji = 1, jpi 
    234                IF ( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
     234               IF( ABS(rab(ji,jj)) < 1._wp ) THEN         ! day duration is less than 24h 
    235235                  rscal(ji,jj) = 0.0_wp 
    236                   IF ( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN      ! day time in one part 
     236                  IF( rdawn_dcy(ji,jj) < rdusk_dcy(ji,jj) ) THEN      ! day time in one part 
    237237                     IF( (rdusk_dcy(ji,jj) - rdawn_dcy(ji,jj) ) .ge. 0.001_wp ) THEN 
    238238                        rscal(ji,jj) = fintegral(rdawn_dcy(ji,jj), rdusk_dcy(ji,jj), raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
     
    247247                  ENDIF 
    248248               ELSE 
    249                   IF ( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day 
     249                  IF( raa(ji,jj) > rbb(ji,jj) ) THEN         ! 24h day 
    250250                     rscal(ji,jj) = fintegral(0._wp, 1._wp, raa(ji,jj), rbb(ji,jj), rcc(ji,jj)) 
    251251                     rscal(ji,jj) = 1._wp / rscal(ji,jj) 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcice_cice.F90

    r11831 r11963  
    132132         IF      ( ksbc == jp_flx ) THEN 
    133133            CALL cice_sbc_force(kt) 
    134          ELSE IF ( ksbc == jp_purecpl ) THEN 
     134         ELSE IF( ksbc == jp_purecpl ) THEN 
    135135            CALL sbc_cpl_ice_flx( fr_i ) 
    136136         ENDIF 
     
    140140         CALL cice_sbc_out ( kt, ksbc ) 
    141141 
    142          IF ( ksbc == jp_purecpl )  CALL cice_sbc_hadgam(kt+1) 
     142         IF( ksbc == jp_purecpl )  CALL cice_sbc_hadgam(kt+1) 
    143143 
    144144      ENDIF                                          ! End sea-ice time step only 
     
    168168      ! there is no restart file. 
    169169      ! Values from a CICE restart file would overwrite this 
    170       IF ( .NOT. ln_rstart ) THEN     
     170      IF( .NOT. ln_rstart ) THEN     
    171171         CALL nemo2cice( tsn(:,:,1,jp_tem) , sst , 'T' , 1.)  
    172172      ENDIF   
     
    177177 
    178178! Do some CICE consistency checks 
    179       IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN 
    180          IF ( calc_strair .OR. calc_Tsfc ) THEN 
     179      IF( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN 
     180         IF( calc_strair .OR. calc_Tsfc ) THEN 
    181181            CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=F and calc_Tsfc=F in ice_in' ) 
    182182         ENDIF 
    183       ELSEIF (ksbc == jp_blk) THEN 
    184          IF ( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN 
     183      ELSEIF(ksbc == jp_blk) THEN 
     184         IF( .NOT. (calc_strair .AND. calc_Tsfc) ) THEN 
    185185            CALL ctl_stop( 'STOP', 'cice_sbc_init : Forcing option requires calc_strair=T and calc_Tsfc=T in ice_in' ) 
    186186         ENDIF 
     
    202202 
    203203      CALL cice2nemo(aice,fr_i, 'T', 1. ) 
    204       IF ( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN 
     204      IF( (ksbc == jp_flx) .OR. (ksbc == jp_purecpl) ) THEN 
    205205         DO jl=1,ncat 
    206206            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. ) 
     
    297297! forced and coupled case  
    298298 
    299       IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN 
     299      IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN 
    300300 
    301301         ztmpn(:,:,:)=0.0 
     
    322322 
    323323! Surface downward latent heat flux (CI_5) 
    324          IF (ksbc == jp_flx) THEN 
     324         IF(ksbc == jp_flx) THEN 
    325325            DO jl=1,ncat 
    326326               ztmpn(:,:,jl)=qla_ice(:,:,1)*a_i(:,:,jl) 
     
    332332            DO jj=1,jpj 
    333333               DO ji=1,jpi 
    334                   IF (fr_i(ji,jj).eq.0.0) THEN 
     334                  IF(fr_i(ji,jj).eq.0.0) THEN 
    335335                     DO jl=1,ncat 
    336336                        ztmpn(ji,jj,jl)=0.0 
     
    351351! GBM conductive flux through ice (CI_6) 
    352352!  Convert to GBM 
    353             IF (ksbc == jp_flx) THEN 
     353            IF(ksbc == jp_flx) THEN 
    354354               ztmp(:,:) = botmelt(:,:,jl)*a_i(:,:,jl) 
    355355            ELSE 
     
    360360! GBM surface heat flux (CI_7) 
    361361!  Convert to GBM 
    362             IF (ksbc == jp_flx) THEN 
     362            IF(ksbc == jp_flx) THEN 
    363363               ztmp(:,:) = (topmelt(:,:,jl)+botmelt(:,:,jl))*a_i(:,:,jl)  
    364364            ELSE 
     
    368368         ENDDO 
    369369 
    370       ELSE IF (ksbc == jp_blk) THEN 
     370      ELSE IF(ksbc == jp_blk) THEN 
    371371 
    372372! Pass bulk forcing fields to CICE (which will calculate heat fluxes etc itself) 
     
    546546! Freshwater fluxes  
    547547 
    548       IF (ksbc == jp_flx) THEN 
     548      IF(ksbc == jp_flx) THEN 
    549549! Note that emp from the forcing files is evap*(1-aice)-(tprecip-aice*sprecip) 
    550550! What we want here is evap*(1-aice)-tprecip*(1-aice) hence manipulation below 
     
    552552! Better to use evap and tprecip? (but for now don't read in evap in this case) 
    553553         emp(:,:)  = emp(:,:)+fr_i(:,:)*(tprecip(:,:)-sprecip(:,:)) 
    554       ELSE IF (ksbc == jp_blk) THEN 
     554      ELSE IF(ksbc == jp_blk) THEN 
    555555         emp(:,:)  = (1.0-fr_i(:,:))*emp(:,:)         
    556       ELSE IF (ksbc == jp_purecpl) THEN 
     556      ELSE IF(ksbc == jp_purecpl) THEN 
    557557! emp_tot is set in sbc_cpl_ice_flx (called from cice_sbc_in above)  
    558558! This is currently as required with the coupling fields from the UM atmosphere 
     
    584584! Scale qsr and qns according to ice fraction (bulk formulae only) 
    585585 
    586       IF (ksbc == jp_blk) THEN 
     586      IF(ksbc == jp_blk) THEN 
    587587         qsr(:,:)=qsr(:,:)*(1.0-fr_i(:,:)) 
    588588         qns(:,:)=qns(:,:)*(1.0-fr_i(:,:)) 
    589589      ENDIF 
    590590! Take into account snow melting except for fully coupled when already in qns_tot 
    591       IF (ksbc == jp_purecpl) THEN 
     591      IF(ksbc == jp_purecpl) THEN 
    592592         qsr(:,:)= qsr_tot(:,:) 
    593593         qns(:,:)= qns_tot(:,:) 
     
    624624 
    625625      CALL cice2nemo(aice,fr_i,'T', 1. ) 
    626       IF ( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN 
     626      IF( (ksbc == jp_flx).OR.(ksbc == jp_purecpl) ) THEN 
    627627         DO jl=1,ncat 
    628628            CALL cice2nemo(aicen(:,:,jl,:),a_i(:,:,jl), 'T', 1. ) 
     
    879879!     B. Gather pn into global array (png) 
    880880 
    881       IF ( jpnij > 1) THEN 
     881      IF( jpnij > 1) THEN 
    882882         CALL mppsync 
    883883         CALL mppgather (pn,0,png)  
     
    892892! (may be OK but not 100% sure) 
    893893 
    894       IF (nproc==0) THEN      
     894      IF(nproc==0) THEN      
    895895!        pcg(:,:)=0.0 
    896896         DO jn=1,jpnij 
     
    10151015! the lbclnk call on pn will replace these with sensible values 
    10161016 
    1017       IF (nproc==0) THEN 
     1017      IF(nproc==0) THEN 
    10181018         png(:,:,:)=0.0 
    10191019         DO jn=1,jpnij 
     
    10281028!     C. Scatter png into NEMO field (pn) for each processor 
    10291029 
    1030       IF ( jpnij > 1) THEN 
     1030      IF( jpnij > 1) THEN 
    10311031         CALL mppsync 
    10321032         CALL mppscatter (png,0,pn)  
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcisf.F90

    r11831 r11963  
    303303      ! 
    304304      ! Allocate public variable 
    305       IF ( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 
     305      IF( sbc_isf_alloc()  /= 0 )         CALL ctl_stop( 'STOP', 'sbc_isf : unable to allocate arrays' ) 
    306306      ! 
    307307      ! initialisation 
     
    440440            !! Initialize arrays to 0 (each step) 
    441441            zt_sum = 0.e0_wp 
    442             IF ( ik > 1 ) THEN 
     442            IF( ik > 1 ) THEN 
    443443               ! 1. -----------the average temperature between 200m and 600m --------------------- 
    444444               DO jk = misfkt(ji,jj),misfkb(ji,jj) 
     
    459459            ELSE 
    460460               qisf(ji,jj) = 0._wp   ;   fwfisf(ji,jj) = 0._wp 
    461             END IF 
     461            ENDIF 
    462462         END DO 
    463463      END DO 
     
    496496      ! coeficient for linearisation of potential tfreez 
    497497      ! Crude approximation for pressure (but commonly used) 
    498       IF ( l_useCT ) THEN   ! linearisation from Jourdain et al. (2017) 
     498      IF( l_useCT ) THEN   ! linearisation from Jourdain et al. (2017) 
    499499         zlamb1 =-0.0564_wp 
    500500         zlamb2 = 0.0773_wp 
     
    558558                  ! compute s freeze 
    559559                  zsfrz=(-zbqe-SQRT(zdis))*zaqer 
    560                   IF ( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 
     560                  IF( zsfrz < 0.0_wp ) zsfrz=(-zbqe+SQRT(zdis))*zaqer 
    561561 
    562562                  ! compute t freeze (eq. 22) 
     
    578578 
    579579         ! define if we need to iterate (nn_gammablk 0/1 do not need iteration) 
    580          IF ( nn_gammablk <  2 ) THEN ; lit = .FALSE. 
     580         IF( nn_gammablk <  2 ) THEN ; lit = .FALSE. 
    581581         ELSE                            
    582582            ! check total number of iteration 
    583             IF (nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
     583            IF(nit >= 100) THEN ; CALL ctl_stop( 'STOP', 'sbc_isf_hol99 : too many iteration ...' ) 
    584584            ELSE                 ; nit = nit + 1 
    585             END IF 
     585            ENDIF 
    586586 
    587587            ! compute error between 2 iterations 
    588588            ! if needed save gammat and compute zhtflx_b for next iteration 
    589589            zerr = MAXVAL(ABS(zhtflx-zhtflx_b)) 
    590             IF ( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 
     590            IF( zerr <= 0.01_wp ) THEN ; lit = .FALSE. 
    591591            ELSE                        ; zhtflx_b(:,:) = zhtflx(:,:) 
    592             END IF 
    593          END IF 
     592            ENDIF 
     593         ENDIF 
    594594      END DO 
    595595      ! 
     
    718718                  pgt(ji,jj) = zustar(ji,jj) / (zgturb + zgmolet) 
    719719                  pgs(ji,jj) = zustar(ji,jj) / (zgturb + zgmoles) 
    720                END IF 
     720               ENDIF 
    721721            END DO 
    722722         END DO 
     
    757757               ! determine the deepest level influenced by the boundary layer 
    758758               DO jk = ikt+1, mbku(ji,jj) 
    759                   IF ( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
     759                  IF( (SUM(e3u_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (umask(ji,jj,jk) == 1) ) ikb = jk 
    760760               END DO 
    761761               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3u_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     
    789789               ! determine the deepest level influenced by the boundary layer 
    790790               DO jk = ikt+1, mbkv(ji,jj) 
    791                   IF ( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
     791                  IF( (SUM(e3v_n(ji,jj,ikt:jk-1)) < zhisf_tbl(ji,jj)) .AND. (vmask(ji,jj,jk) == 1) ) ikb = jk 
    792792               END DO 
    793793               zhisf_tbl(ji,jj) = MIN(zhisf_tbl(ji,jj), SUM(e3v_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     
    869869               ! determine the deepest level influenced by the boundary layer 
    870870               DO jk = ikt, mbkt(ji,jj) 
    871                   IF ( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
     871                  IF( (SUM(e3t_n(ji,jj,ikt:jk-1)) .LT. rhisf_tbl(ji,jj)) .AND. (tmask(ji,jj,jk) == 1) ) ikb = jk 
    872872               END DO 
    873873               rhisf_tbl(ji,jj) = MIN(rhisf_tbl(ji,jj), SUM(e3t_n(ji,jj,ikt:ikb)))  ! limit the tbl to water thickness. 
     
    879879            END DO 
    880880         END DO 
    881       END IF  
     881      ENDIF  
    882882      ! 
    883883      !==   ice shelf melting distributed over several levels   ==! 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcmod.F90

    r11831 r11963  
    243243         fwfisf  (:,:)   = 0._wp   ;   risf_tsc  (:,:,:) = 0._wp 
    244244         fwfisf_b(:,:)   = 0._wp   ;   risf_tsc_b(:,:,:) = 0._wp 
    245       END IF 
     245      ENDIF 
    246246      IF( nn_ice == 0 ) THEN        !* No sea-ice in the domain : ice fraction is always zero 
    247247         IF( nn_components /= jp_iam_opa )   fr_i(:,:) = 0._wp    ! except for OPA in SAS-OPA coupled case 
     
    399399         emp_b (:,:) = emp (:,:) 
    400400         sfx_b (:,:) = sfx (:,:) 
    401          IF ( ln_rnf ) THEN 
     401         IF( ln_rnf ) THEN 
    402402            rnf_b    (:,:  ) = rnf    (:,:  ) 
    403403            rnf_tsc_b(:,:,:) = rnf_tsc(:,:,:) 
     
    437437      IF( ln_mixcpl )          CALL sbc_cpl_rcv   ( kt, nn_fsbc, nn_ice )   ! forced-coupled mixed formulation after forcing 
    438438      ! 
    439       IF ( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )      ! Wind stress provided by waves  
     439      IF( ln_wave .AND. (ln_tauwoc .OR. ln_tauw) ) CALL sbc_wstress( )      ! Wind stress provided by waves  
    440440      ! 
    441441      !                                            !==  Misc. Options  ==! 
     
    471471!!$!clem: it looks like it is necessary for the north fold (in certain circumstances). Don't know why. 
    472472!!$      CALL lbc_lnk( 'sbcmod', emp, 'T', 1. ) 
    473       IF ( ll_wd ) THEN     ! If near WAD point limit the flux for now 
     473      IF( ll_wd ) THEN     ! If near WAD point limit the flux for now 
    474474         zthscl = atanh(rn_wd_sbcfra)                     ! taper frac default is .999  
    475475         zwdht(:,:) = sshn(:,:) + ht_0(:,:) - rn_wdmin1   ! do this calc of water 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbcrnf.F90

    r11831 r11963  
    439439         !                                      !    - mixed upstream-centered (ln_traadv_cen2=T) 
    440440         ! 
    441          IF ( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   & 
     441         IF( ln_rnf_depth )   CALL ctl_warn( 'sbc_rnf_init: increased mixing turned on but effects may already',   & 
    442442            &                                              'be spread through depth by ln_rnf_depth'               ) 
    443443         ! 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/sbctide.F90

    r10068 r11963  
    7272         ! Temporarily set nsec_day to beginning of day. 
    7373         nsec_day_orig = nsec_day 
    74          IF ( nsec_day /= NINT(0.5_wp * rdt) ) THEN  
     74         IF( nsec_day /= NINT(0.5_wp * rdt) ) THEN  
    7575            kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt 
    7676            nsec_day = NINT(0.5_wp * rdt) 
  • NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/src/OCE/SBC/tideini.F90

    r11831 r11963  
    6868      ! 
    6969      IF( ln_tide ) THEN 
    70          IF (lwp) THEN 
     70         IF(lwp) THEN 
    7171            WRITE(numout,*) 
    7272            WRITE(numout,*) 'tide_init : Initialization of the tidal components' 
     
    127127      kt_tide = nit000 
    128128      ! 
    129       IF (.NOT.ln_scal_load ) rn_scal_load = 0._wp 
     129      IF(.NOT.ln_scal_load ) rn_scal_load = 0._wp 
    130130      ! 
    131131   END SUBROUTINE tide_init 
Note: See TracChangeset for help on using the changeset viewer.