C $Header: /home/ubuntu/mnt/e9_copy/MITgcm_contrib/AITCZ/code/mitphys_do_inphys.F,v 1.1 2003/08/20 15:24:59 czaja Exp $ C $Name: $ #include "CPP_OPTIONS.h" SUBROUTINE MITPHYS_INIT( myThid ) C /==================================================================\ C | S/R MITPHYS_INIT | C |==================================================================| C | Interface between MITPHYS atmospheric physics package and the | C | dynamical model for initialisation. | C \==================================================================/ IMPLICIT NONE C -------------- Global variables ------------------------------------ #include "atparam.h" #include "atparam1.h" #include "EEPARAMS.h" #include "PARAMS.h" #include "DYNVARS.h" #include "GRID.h" #include "MITPHYS_DIAGS.h" #include "MITPHYS_PARAMS.h" INTEGER NGP INTEGER NLON INTEGER NLAT INTEGER NLEV PARAMETER ( NLON=IX, NLAT=IL, NLEV=KX, NGP=NLON*NLAT ) #include "com_mitphysvar.h" #include "com_forcing1.h" C == Routine arguments == C myThid - Number of this instance INTEGER myThid #ifdef ALLOW_MITPHYS C == Local variables == C FSG - Cell mid-point in vertical C HSG - Cell face in vertical C RLAT - Call mid-point latitude C pGround - Lower boundary pressure C J, K, bi,bj - Loop counters _RL FSG( Nr) _RL HSG( Nr+1) _RL RLAT(sNy) _RL pGround INTEGER I,I2,J, K INTEGER Katm INTEGER bi,bj _RL LAT_CYCLE, DIST _RS FMASK2(1-OLx:sNx+OLx,1-OLy:sNy+OLy,nSx,nSy) pground = 1012.5D2 c omp: Single thread per processor: bi = 1 bj = 1 DO K=1,Nr Katm = K FSG(Katm ) = rC(K)/pGround HSG(Katm+1) = rF(K)/pGround ENDDO HSG(1) = rF(Nr+1)/pGround DO J = 1, sNy DO I = 1, sNx I2 = (sNx)*(J-1)+I if (UsingSphericalPolarGrid) THEN C latitude and logitude, assuming a spherical coordinate system LAT_G(I2) = yC(I,J,bi,bj) LON_G(I2) = xC(I,J,bi,bj) ELSE C For Cartesian Grid, need to specify - equatorial value for now... LAT_G(I2) = 0.0 LON_G(I2) = 0.0 END IF CBMFG1 (I2)= 0.0 END DO END DO C DO J=1,sNy C RLAT(J) = yC(1,J,1,1)*deg2rad C ENDDO CALL MITPHYS_READPARMS(myThid) C 2. SST / STL C modif omp: the background SST is defined thourgh the MITPHYS namelist C modif acz: introduce analytical STL C omp: for spherical geometry IF (UsingSphericalPolarGrid) THEN LAT_CYCLE = 1.0 DO J=1,sNy DO I=1,sNx I2 = (sNx)*(J-1)+I SST_REF(I2) = SST_BACK STL_REF(I2) = SST_BACK C (acz): NTA positive SST anomaly dist = SQRT ( ( LAT_G(I2) - LAT_PERT_NTA) ** 2) if (dist .LE. DEL_LAT_NTA) then SST_REF(I2) = SST_REF(I2) + DELT_PERT_NTA & * cos ( pi * dist / 2./DEL_LAT_NTA ) **2 end if C Equator to pole SST gradient SST_REF(I2) = SST_REF(I2) - DELT_EQ_PO * & sin( LAT_G(I2) * pi / 180. ) **2 C Equator to northern boundary SST gradient SST_REF(I2) = SST_REF(I2) - DELT_EQ_BND * & sin( LAT_G(I2) / abs(Phimin)) **2 C Lindzen and Hou type gradient: SST_REF(I2) = SST_REF(I2) - DELT_LH * ( & 3 * ( sin ( PI * LAT_G(I2) / 180.d0 ) & - sin(PI * LAT_CYCLE * LAT_LH / 180.d0) & ) ** 2) CC (acz) The SST used for TOSA1: CC Zonally symmetric SST perturbation dist = SQRT ( ( LAT_G(I2) - LAT_PERT*LAT_CYCLE) ** 2) if (dist .LE. DEL_LAT) then SST_REF(I2) = SST_REF(I2) + DELT_PERT & * cos ( pi * dist / 2./DEL_LAT ) **2 end if CC (acz) STL used for TOSA1: CC CC 1 - Warm ring centered at the equator with same CC meridional extent (DEL_LAT) than the above (TOSA1) SST CC The profile is further time-dependent: C strongest at equinox, zero at solstice C so it is set at zero for this initialization dist = ABS ( LAT_G(I2) ) if (dist .LE. DEL_LAT) then STL_REF(I2) = STL_REF(I2) + DELT_STL_EQU & * cos ( pi * dist / 2./DEL_LAT ) **2 & * cos ( 2 * pi * (0. / 24./ 3600.) & / 360. + 0.5 * pi ) **2 end if C 2- Opposite sign variations north and south of equator C North is warm at t=0, dist = ABS ( LAT_G(I2) ) if (dist .LE. DEL_LAT_STL) then STL_REF(I2) = STL_REF(I2) + DELT_STL_HEM & * sin ( pi * LAT_G(I2) / DEL_LAT_STL ) end if ENDDO ENDDO ELSE C Cartesian Grid CC (acz) nothing for the moment ... END IF C 3. SOIL MOISTURE C (acz): start with full bucket C for spherical geometry IF (UsingSphericalPolarGrid) THEN DO J=1,sNy DO I=1,sNx I2 = (sNx)*(J-1)+I BUCKET(I2) = BKT0 END DO END DO CC(acz) CC CALL WRITE_FLD_XY_RS('BKTini','0000000000',BUCKET,0,myThid) CC write(6,*) 'end of do_inphys: BKT0=', BKT0 ELSE C Cartesian Grid: nothing is done for now... END IF C 4. LAND-SEA MASK C (acz): FMASK1 is defined in com_forcing1.h It is determined analytically here C FMASK1 = 0 over ocean C FMASK1 = 1 over land C for spherical geometry IF (UsingSphericalPolarGrid) THEN DO J=1,sNy DO I=1,sNx I2 = (sNx)*(J-1)+I FMASK1(I2) = 0 cc ATL1 experiment: AMERICA - ATLANTIC - AFRICA cc each block is defined within two meridians cc bounded by 30S and 60N (maximum northward extent cc of the model) FMASK1(I2) = 0.0 if ( (LON_G(I2) .GE. 105.) .AND. (LON_G(I2) .LE. 150.) .AND. (LAT_G(I2) .GE. -30) ) then FMASK1(I2) = 1.0 end if if ( (LON_G(I2) .GE. 210.) .AND. (LON_G(I2) .LE. 255.) .AND. (LAT_G(I2) .GE. -30) ) then FMASK1(I2) = 1.0 end if END DO END DO ELSE C Cartesian Grid: nothing is done for now... END IF C 5. QFLUX C C for spherical geometry IF (UsingSphericalPolarGrid) THEN DO J=1,sNy DO I=1,sNx I2 = (sNx)*(J-1)+I cc QSTAR: `equinox', centered on the equator cc with meridional scale set by DELL_STAR and strength cc by DELQ_STAR. This is added to a background cst value QSTAR_BACK QSTAR(I2) = QSTAR_BACK dist = ABS ( LAT_G(I2) ) if (dist .LE. DELL_STAR) then QSTAR(I2) = QSTAR(I2) + DELQ_STAR * cos ( pi * dist / 2./DELL_STAR ) **2 end if cc cc QOCEAN: MOC and / or ENSO forcing cc if (FMASK1(I2) .EQ. 1) then QOCEAN(I2) = 0. else QOCEAN(I2) = 0. cc MOC: meridional dipole across the equator only over the Atlantic if ( (LON_G(I2) .GE. 150.) .AND. (LON_G(I2) .LE. 210.) ) then dist = ABS ( LAT_G(I2) ) if (dist .LE. DELL_OCEAN) then QOCEAN(I2) = DELQ_OCEAN * sin ( pi * LAT_G(I2) / DELL_OCEAN ) else QOCEAN(I2) = 0. end if end if cc ENSO: Nino3 monopole in the `tropical Pacific' cc with meridional scale DELL_OCEAN, non zero between longitudes cc LOE_OCEAN and LOW_OCEAN, and amplitude DELQ_OCEAN cc cc !!!!!! check that QOCEAN is not reset when entering here... cc + define a new DELQ_OCEAN for ENSO... cc cc if ( (LON_G(I2) .LE. LOW_OCEAN) .AND. (LON_G(I2) .GE. LOE_OCEAN) ) then cc QOCEAN(I2) = DELQ_OCEAN * sin( pi*(LON_G(I2)-LOE_OCEAN) cc & / (LOW_OCEAN-LOE_OCEAN) ) cc end if cc dist = ABS ( LAT_G(I2) ) cc if (dist .LE. DELL_OCEAN) then cc QOCEAN(I2) = QOCEAN(I2) * cos ( pi * dist / 2./DELL_OCEAN ) **2 cc else cc QOCEAN(I2) = 0. cc end if end if QFLUX(I2) = QSTAR(I2) + QOCEAN(I2) END DO END DO CC(acz) CC CALL WRITE_FLD_XY_RS('QFLUX','0000000000',QFLUX,0,myThid) ELSE C Cartesian Grid: nothing is done for now... END IF C DO J=1,sNy C DO I=1,sNx C I2 = (sNx)*(J-1)+I C SST1 (I2)= SST_REF(I2) C CBMF_DYN(I,J,bi,bj) = CBMFG1(I2) C SST_DYN(I,J,bi,bj) = SST_REF(I2) C END DO C END DO C C _EXCH_XY_R8(CBMF_DYN,myThid) C _EXCH_XY_R8(SST_DYN,myThid) c -- for the MIT physics, this is not necessary anymore: ccc CALL INPHYS( FSG, HSG, RLAT ) c -- #ifdef ALLOW_TIMEAVE C Initialise diagnostic counters C (these are cleared on model starti i.e. not C loaded from history file for now ). DO bj = myByLo(myThid), myByHi(myThid) DO bi = myBxLo(myThid), myBxHi(myThid) CALL TIMEAVE_RESET(TSWtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(TLWtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(SSWtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(SLWtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(SINStave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(BKTtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(SHFtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(LHFtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(PRECNVtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(PRECLStave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(CLOUDCtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(SSTtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(PRWtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(TSW0tave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(TLW0tave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(CBMFtave, 1, bi, bj, myThid) CALL TIMEAVE_RESET(RHtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(CLDFtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(CLDQtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(CLDQCtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(RCtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(CRLWtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(CRSWtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(MUPtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(MDNtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(MDN0tave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(ENTtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(DETtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(Utave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(Vtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(Wtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(Ttave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(Qtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(PHItave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(SXtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(HXtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(DTCNVtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(DTLSCtave, Nr, bi, bj, myThid) CALL TIMEAVE_RESET(DTPBLtave, Nr, bi, bj, myThid) DO K = 1, Nr MITPHYS_TimeAve(K,bi,bj) = 0. ENDDO ENDDO ENDDO #endif /* ALLOW_TIMEAVE */ #endif /* ALLOW_MITPHYS */ RETURN END