Ignore:
Timestamp:
2017-08-09T14:56:31+02:00 (3 years ago)
Author:
dford
Message:

Add logchl assimilation, with nitrogen balancing for HadOCC.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90

    r8400 r8428  
    2222   !!   ssh_asm_inc    : Apply the SSH increment 
    2323   !!   seaice_asm_inc : Apply the seaice increment 
     24   !!   logchl_asm_inc : Apply the logchl increment 
    2425   !!---------------------------------------------------------------------- 
    2526   USE wrk_nemo         ! Memory Allocation 
     
    4445#endif 
    4546   USE sbc_oce          ! Surface boundary condition variables. 
     47   USE zdfmxl, ONLY :  &   
     48!   &  hmld_tref,       &    
     49#if defined key_karaml 
     50   &  hmld_kara,       & 
     51   &  ln_kara,         & 
     52#endif    
     53   &  hmld,            &  
     54   &  hmlp,            & 
     55   &  hmlpt  
     56#if defined key_top 
     57   USE trc, ONLY: & 
     58      & trn,      & 
     59      & trb 
     60   USE par_trc, ONLY: & 
     61      & jptra 
     62#endif 
     63#if defined key_medusa && defined key_foam_medusa 
     64   USE asmlogchlbal_medusa, ONLY: & 
     65      & asm_logchl_bal_medusa 
     66   USE par_medusa 
     67#elif defined key_hadocc 
     68   USE asmlogchlbal_hadocc, ONLY: & 
     69      & asm_logchl_bal_hadocc 
     70   USE par_hadocc 
     71#endif 
    4672 
    4773   IMPLICIT NONE 
     
    5480   PUBLIC   ssh_asm_inc    !: Apply the SSH increment 
    5581   PUBLIC   seaice_asm_inc !: Apply the seaice increment 
     82   PUBLIC   logchl_asm_inc !: Apply the logchl increment 
    5683 
    5784#if defined key_asminc 
     
    6188#endif 
    6289   LOGICAL, PUBLIC :: ln_bkgwri = .FALSE.      !: No output of the background state fields 
     90   LOGICAL, PUBLIC :: ln_balwri = .FALSE.      !: No output of the assimilation balancing increments 
    6391   LOGICAL, PUBLIC :: ln_asmiau = .FALSE.      !: No applying forcing with an assimilation increment 
    6492   LOGICAL, PUBLIC :: ln_asmdin = .FALSE.      !: No direct initialization 
     
    6795   LOGICAL, PUBLIC :: ln_sshinc = .FALSE.      !: No sea surface height assimilation increment 
    6896   LOGICAL, PUBLIC :: ln_seaiceinc             !: No sea ice concentration increment 
     97   LOGICAL, PUBLIC :: ln_logchlinc = .FALSE.   !: No log10(chlorophyll) increment 
     98   LOGICAL, PUBLIC :: ln_logchlbal = .FALSE.   !: Don't apply multivariate log10(chlorophyll) balancing 
    6999   LOGICAL, PUBLIC :: ln_salfix = .FALSE.      !: Apply minimum salinity check 
    70100   LOGICAL, PUBLIC :: ln_temnofreeze = .FALSE. !: Don't allow the temperature to drop below freezing 
     
    91121   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   ssh_bkg, ssh_bkginc   ! Background sea surface height and its increment 
    92122   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   seaice_bkginc         ! Increment to the background sea ice conc 
     123 
     124   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: logchl_bkginc !: Increment to background logchl 
     125#if defined key_top 
     126   REAL(wp), PUBLIC, DIMENSION(:,:,:,:), ALLOCATABLE :: logchl_balinc  !: Increment to BGC variables from logchl assim 
     127#endif 
     128   REAL(wp) :: rn_maxchlinc = -999.0  !: maximum absolute non-log chlorophyll increment from logchl assimilation 
     129                                      !: <= 0 implies no maximum applied (switch turned off) 
     130                                      !:  > 0 implies maximum absolute chl increment capped at this value 
     131 
     132   INTEGER :: mld_choice_bgc    = 5   !: choice of mld criteria to use for biogeochemistry assimilation 
     133                                      !: 1) hmld      - Turbocline/mixing depth                           [W points] 
     134                                      !: 2) hmlp      - Density criterion (0.01 kg/m^3 change from 10m)   [W points] 
     135                                      !: 3) hmld_kara - Kara MLD                                          [Interpolated] 
     136                                      !: 4) hmld_tref - Temperature criterion (0.2 K change from surface) [T points] 
     137                                      !: 5) hmlpt     - Density criterion (0.01 kg/m^3 change from 10m)   [T points] 
    93138 
    94139   !! * Substitutions 
     
    133178      REAL(wp), POINTER, DIMENSION(:,:) ::   hdiv   ! 2D workspace 
    134179      !! 
    135       NAMELIST/nam_asminc/ ln_bkgwri,                                      & 
     180      NAMELIST/nam_asminc/ ln_bkgwri, ln_balwri,                           & 
    136181         &                 ln_trainc, ln_dyninc, ln_sshinc,                & 
     182         &                 ln_logchlinc, ln_logchlbal,                     & 
    137183         &                 ln_asmdin, ln_asmiau,                           & 
    138184         &                 nitbkg, nitdin, nitiaustr, nitiaufin, niaufn,   & 
    139185         &                 ln_salfix, salfixmin, nn_divdmp,                & 
    140          &                 ln_seaiceinc, ln_temnofreeze 
     186         &                 ln_seaiceinc, ln_temnofreeze,                   & 
     187         &                 mld_choice_bgc, rn_maxchlinc 
    141188      !!---------------------------------------------------------------------- 
    142189 
     
    163210         WRITE(numout,*) '   Namelist namasm : set assimilation increment parameters' 
    164211         WRITE(numout,*) '      Logical switch for writing out background state          ln_bkgwri = ', ln_bkgwri 
     212         WRITE(numout,*) '      Logical switch for writing out balancing increments      ln_balwri = ', ln_balwri 
    165213         WRITE(numout,*) '      Logical switch for applying tracer increments            ln_trainc = ', ln_trainc 
    166214         WRITE(numout,*) '      Logical switch for applying velocity increments          ln_dyninc = ', ln_dyninc 
     
    168216         WRITE(numout,*) '      Logical switch for Direct Initialization (DI)            ln_asmdin = ', ln_asmdin 
    169217         WRITE(numout,*) '      Logical switch for applying sea ice increments        ln_seaiceinc = ', ln_seaiceinc 
     218         WRITE(numout,*) '      Logical switch for applying logchl increments         ln_logchlinc = ', ln_logchlinc 
     219         WRITE(numout,*) '      Logical switch for logchl multivariate balancing      ln_logchlbal = ', ln_logchlbal 
    170220         WRITE(numout,*) '      Logical switch for Incremental Analysis Updating (IAU)   ln_asmiau = ', ln_asmiau 
    171221         WRITE(numout,*) '      Timestep of background in [0,nitend-nit000-1]            nitbkg    = ', nitbkg 
     
    176226         WRITE(numout,*) '      Logical switch for ensuring that the sa > salfixmin      ln_salfix = ', ln_salfix 
    177227         WRITE(numout,*) '      Minimum salinity after applying the increments           salfixmin = ', salfixmin 
     228         WRITE(numout,*) '      Choice of MLD for BGC assimilation                  mld_choice_bgc = ', mld_choice_bgc 
     229         WRITE(numout,*) '      Maximum absolute chlorophyll increment (<=0 = off)    rn_maxchlinc = ', rn_maxchlinc 
    178230      ENDIF 
    179231 
     
    223275 
    224276      IF (      ( ( .NOT. ln_asmdin ).AND.( .NOT. ln_asmiau ) ) & 
    225            .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ) .OR. ( ln_seaiceinc) )) & 
    226          & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc is set to .true.', & 
     277         & .AND.( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ).OR. & 
     278         &        ( ln_logchlinc ) )) & 
     279         & CALL ctl_stop( ' One or more of ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 
     280         &                ' and ln_logchlinc is set to .true.', & 
    227281         &                ' but ln_asmdin and ln_asmiau are both set to .false. :', & 
    228282         &                ' Inconsistent options') 
     
    233287 
    234288      IF ( ( .NOT. ln_trainc ).AND.( .NOT. ln_dyninc ).AND.( .NOT. ln_sshinc ).AND.( .NOT. ln_seaiceinc ) & 
    235          &                     )  & 
    236          & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc and ln_seaiceinc are set to .false. :', & 
     289         & .AND.( .NOT. ln_logchlinc ) )  & 
     290         & CALL ctl_warn( ' ln_trainc, ln_dyninc, ln_sshinc, ln_seaiceinc,', & 
     291         &                ' and ln_logchlinc are set to .false. :', & 
    237292         &                ' The assimilation increments are not applied') 
    238293 
     
    254309         &                ' Background time step for Direct Initialization is outside', & 
    255310         &                ' the cycle interval') 
     311 
     312      IF ( ( ln_balwri ).AND.( .NOT. ( ln_logchlinc ) ) ) THEN 
     313         CALL ctl_warn( ' Balancing increments are only calculated for logchl', & 
     314            &           ' Not assimilating logchl, so ln_balwri will be set to .false.') 
     315         ln_balwri = .FALSE. 
     316      ENDIF 
     317 
     318      IF ( ( ln_logchlbal ).AND.( .NOT. ( ln_logchlinc ) ) ) THEN 
     319         CALL ctl_warn( ' Cannot calculate logchl balancing increments if not assimilating logchl', & 
     320            &           ' Not assimilating logchl, so ln_logchlbal will be set to .false.') 
     321         ln_logchlbal = .FALSE. 
     322      ENDIF 
    256323 
    257324      IF ( nstop > 0 ) RETURN       ! if there are any errors then go no further 
     
    350417      ssh_iau(:,:)    = 0.0 
    351418#endif 
    352       IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) ) THEN 
     419      IF ( ln_logchlinc ) THEN 
     420         ALLOCATE( logchl_bkginc(jpi,jpj) ) 
     421         logchl_bkginc(:,:) = 0.0 
     422#if defined key_top 
     423         ALLOCATE( logchl_balinc(jpi,jpj,jpk,jptra) ) 
     424         logchl_balinc(:,:,:,:) = 0.0 
     425#endif 
     426      ENDIF 
     427      IF ( ( ln_trainc ).OR.( ln_dyninc ).OR.( ln_sshinc ).OR.( ln_seaiceinc ) & 
     428         &  .OR.( ln_logchlinc ) ) THEN 
    353429 
    354430         !-------------------------------------------------------------------- 
     
    422498            ! to allow for differences in masks 
    423499            WHERE( ABS( seaice_bkginc(:,:) ) > 1.0e+10 ) seaice_bkginc(:,:) = 0.0 
     500         ENDIF 
     501 
     502         IF ( ln_logchlinc ) THEN 
     503            CALL iom_get( inum, jpdom_autoglo, 'bckinlogchl', logchl_bkginc(:,:), 1 ) 
     504            ! Apply the masks 
     505            logchl_bkginc(:,:) = logchl_bkginc(:,:) * tmask(:,:,1) 
     506            ! Set missing increments to 0.0 rather than 1e+20 
     507            ! to allow for differences in masks 
     508            WHERE( ABS( logchl_bkginc(:,:) ) > 1.0e+10 ) logchl_bkginc(:,:) = 0.0 
    424509         ENDIF 
    425510 
     
    11551240 
    11561241   END SUBROUTINE seaice_asm_inc 
     1242 
     1243   SUBROUTINE logchl_asm_inc( kt ) 
     1244      !!---------------------------------------------------------------------- 
     1245      !!                    ***  ROUTINE logchl_asm_inc  *** 
     1246      !!           
     1247      !! ** Purpose : Apply the chlorophyll assimilation increments. 
     1248      !! 
     1249      !! ** Method  : Calculate increments to state variables using nitrogen 
     1250      !!              balancing. 
     1251      !!              Direct initialization or Incremental Analysis Updating. 
     1252      !! 
     1253      !! ** Action  :  
     1254      !!---------------------------------------------------------------------- 
     1255      INTEGER, INTENT(IN) :: kt   ! Current time step 
     1256      ! 
     1257      INTEGER  :: jk              ! Loop counter 
     1258      INTEGER  :: it              ! Index 
     1259      REAL(wp) :: zincwgt         ! IAU weight for current time step 
     1260      REAL(wp) :: zincper         ! IAU interval in seconds 
     1261      !!---------------------------------------------------------------------- 
     1262 
     1263      IF ( kt <= nit000 ) THEN 
     1264 
     1265         zincper = (nitiaufin_r - nitiaustr_r + 1) * rdt 
     1266 
     1267#if defined key_medusa && defined key_foam_medusa 
     1268         !CALL asm_logchl_bal_medusa() 
     1269         CALL ctl_stop( 'Attempting to assimilate logchl into MEDUSA, ', & 
     1270            &           'but not fully implemented yet' ) 
     1271#elif defined key_hadocc 
     1272         CALL asm_logchl_bal_hadocc( logchl_bkginc, zincper, mld_choice_bgc, & 
     1273            &                        rn_maxchlinc, ln_logchlbal, logchl_balinc ) 
     1274#else 
     1275         CALL ctl_stop( 'Attempting to assimilate logchl, ', & 
     1276            &           'but not defined a biogeochemical model' ) 
     1277#endif 
     1278 
     1279      ENDIF 
     1280 
     1281      IF ( ln_asmiau ) THEN 
     1282 
     1283         !-------------------------------------------------------------------- 
     1284         ! Incremental Analysis Updating 
     1285         !-------------------------------------------------------------------- 
     1286 
     1287         IF ( ( kt >= nitiaustr_r ).AND.( kt <= nitiaufin_r ) ) THEN 
     1288 
     1289            it = kt - nit000 + 1 
     1290            zincwgt = wgtiau(it)      ! IAU weight for the current time step 
     1291            ! note this is not a tendency so should not be divided by rdt 
     1292 
     1293            IF(lwp) THEN 
     1294               WRITE(numout,*)  
     1295               WRITE(numout,*) 'logchl_asm_inc : logchl IAU at time step = ', & 
     1296                  &  kt,' with IAU weight = ', wgtiau(it) 
     1297               WRITE(numout,*) '~~~~~~~~~~~~' 
     1298            ENDIF 
     1299 
     1300            ! Update the biogeochemical variables 
     1301            ! Add directly to trn and trb, rather than to tra, as not a tendency 
     1302#if defined key_medusa && defined key_foam_medusa 
     1303            DO jk = 1, jpkm1 
     1304               trn(:,:,jk,jp_msa0:jp_msa1) = trn(:,:,jk,jp_msa0:jp_msa1) + & 
     1305                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
     1306               trb(:,:,jk,jp_msa0:jp_msa1) = trb(:,:,jk,jp_msa0:jp_msa1) + & 
     1307                  &                          logchl_balinc(:,:,jk,jp_msa0:jp_msa1) * zincwgt 
     1308            END DO 
     1309#elif defined key_hadocc 
     1310            DO jk = 1, jpkm1 
     1311               trn(:,:,jk,jp_had0:jp_had1) = trn(:,:,jk,jp_had0:jp_had1) + & 
     1312                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 
     1313               trb(:,:,jk,jp_had0:jp_had1) = trb(:,:,jk,jp_had0:jp_had1) + & 
     1314                  &                          logchl_balinc(:,:,jk,jp_had0:jp_had1) * zincwgt 
     1315            END DO 
     1316#endif 
     1317 
     1318            ! Do not deallocate arrays - needed by asm_bal_wri 
     1319            ! which is called at end of model run 
     1320         ENDIF 
     1321 
     1322      ELSEIF ( ln_asmdin ) THEN  
     1323 
     1324         !-------------------------------------------------------------------- 
     1325         ! Direct Initialization 
     1326         !-------------------------------------------------------------------- 
     1327          
     1328         IF ( kt == nitdin_r ) THEN 
     1329 
     1330            neuler = 0                    ! Force Euler forward step 
     1331 
     1332#if defined key_medusa && defined key_foam_medusa 
     1333            ! Initialize the now fields with the background + increment 
     1334            ! Background currently is what the model is initialised with 
     1335            CALL ctl_warn( ' Doing direct initialisation of MEDUSA with chlorophyll assimilation', & 
     1336               &           ' Background state is taken from model rather than background file' ) 
     1337            trn(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) + & 
     1338               &                         logchl_balinc(:,:,:,jp_msa0:jp_msa1) 
     1339            trb(:,:,:,jp_msa0:jp_msa1) = trn(:,:,:,jp_msa0:jp_msa1) 
     1340#elif defined key_hadocc 
     1341            ! Initialize the now fields with the background + increment 
     1342            ! Background currently is what the model is initialised with 
     1343            CALL ctl_warn( ' Doing direct initialisation of HadOCC with chlorophyll assimilation', & 
     1344               &           ' Background state is taken from model rather than background file' ) 
     1345            trn(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) + & 
     1346               &                         logchl_balinc(:,:,:,jp_had0:jp_had1) 
     1347            trb(:,:,:,jp_had0:jp_had1) = trn(:,:,:,jp_had0:jp_had1) 
     1348#endif 
     1349  
     1350            ! Do not deallocate arrays - needed by asm_bal_wri 
     1351            ! which is called at end of model run 
     1352         ENDIF 
     1353         ! 
     1354      ENDIF 
     1355      ! 
     1356   END SUBROUTINE logchl_asm_inc 
    11571357    
    11581358   !!====================================================================== 
Note: See TracChangeset for help on using the changeset viewer.