Changeset 226 for trunk/LMDZ.MARS/libf
- Timestamp:
- Jul 15, 2011, 2:55:17 PM (13 years ago)
- Location:
- trunk/LMDZ.MARS/libf/phymars
- Files:
-
- 10 added
- 1 deleted
- 6 edited
- 1 moved
-
inifis.F (modified) (5 diffs)
-
meso_callkeys.h (deleted)
-
meso_inc (added)
-
meso_inc/meso_dimphys.h_ref (moved) (moved from trunk/LMDZ.MARS/libf/phymars/meso_dimphys.h_ref)
-
meso_inc/meso_inc_caps.F (added)
-
meso_inc/meso_inc_ini.F (added)
-
meso_inc/meso_inc_inifisini.F (added)
-
meso_inc/meso_inc_inifisinvar.F (added)
-
meso_inc/meso_inc_inifisvar.F (added)
-
meso_inc/meso_inc_invar.F (added)
-
meso_inc/meso_inc_les.F (added)
-
meso_inc/meso_inc_slope.F (added)
-
meso_inc/meso_inc_var.F (added)
-
meso_inifis.F (modified) (10 diffs)
-
meso_param_slope.F90 (modified) (4 diffs)
-
meso_physiq.F (modified) (36 diffs)
-
physiq.F (modified) (4 diffs)
-
surfdat.h (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/LMDZ.MARS/libf/phymars/inifis.F
r161 r226 1 SUBROUTINE inifis(ngrid,nlayer, 2 $ day_ini,pdaysec,ptimestep, 3 $ plat,plon,parea, 4 $ prad,pg,pr,pcpp) 1 SUBROUTINE inifis( 2 $ ngrid,nlayer 3 $ ,day_ini,pdaysec,ptimestep 4 $ ,plat,plon,parea 5 $ ,prad,pg,pr,pcpp 6 $ ) 5 7 ! 6 8 !======================================================================= … … 55 57 #include "datafile.h" 56 58 57 REAL prad,pg,pr,pcpp,pdaysec,ptimestep 58 59 REAL prad,pg,pr,pcpp,pdaysec 60 61 REAL ptimestep 62 INTEGER day_ini 63 59 64 INTEGER ngrid,nlayer 60 65 REAL plat(ngrid),plon(ngrid),parea(ngridmx) 61 integer day_ini62 66 INTEGER ig,ierr 63 67 … … 75 79 76 80 rad=prad 77 daysec=pdaysec78 dtphys=ptimestep79 81 cpp=pcpp 80 82 g=pg 81 83 r=pr 82 84 rcp=r/cpp 85 daysec=pdaysec 86 dtphys=ptimestep 83 87 84 88 ! -------------------------------------------------------- … … 578 582 ! ------------------------ 579 583 584 ! in 'comgeomfi.h' 580 585 CALL SCOPY(ngrid,plon,1,long,1) 581 586 CALL SCOPY(ngrid,plat,1,lati,1) … … 583 588 totarea=SSUM(ngridmx,area,1) 584 589 590 ! in 'comdiurn.h' 585 591 DO ig=1,ngrid 586 592 sinlat(ig)=sin(plat(ig)) -
trunk/LMDZ.MARS/libf/phymars/meso_inifis.F
r162 r226 1 SUBROUTINE meso_inifis(ngrid,nlayer,nq, 2 $ wdt, 3 $ wday_ini,wdaysec, 4 $ wappel_phys,wecri_phys, 5 $ plat,plon,parea, 6 $ prad,pg,pr,pcpp, 7 $ womeg,wmugaz, 8 $ wyear_day,wperiheli,waphelie,wperi_day,wobliquit, 9 $ wz0,wemin_turb,wlmixmin, 10 $ wemissiv,wemissiceN,wemissiceS,walbediceN,walbediceS, 11 $ wiceradiusN,wiceradiusS,wdtemisiceN,wdtemisiceS, 12 $ walbedodat,wphisfi,wvolcapa, 13 $ wzmea,wzstd,wzsig,wzgam,wzthe, 14 $ wtheta,wpsi) 15 16 c======================================================================= 17 c 18 c CAREFUL: THIS IS A VERSION TO BE USED WITH WRF !!! 19 c 20 c ... CHECK THE ****WRF lines 21 c 22 c======================================================================= 1 SUBROUTINE meso_inifis( 2 $ ngrid,nlayer 3 #ifdef MESOSCALE 4 $ ,nq,wdt,wday_ini,wdaysec,wappel_phys 5 #else 6 $ ,day_ini,pdaysec,ptimestep 7 #endif 8 $ ,plat,plon,parea 9 $ ,prad,pg,pr,pcpp 10 #ifdef MESOSCALE 11 #include "meso_inc/meso_inc_inifisinvar.F" 12 #endif 13 $ ) 23 14 ! 24 15 !======================================================================= … … 59 50 ! ------------- 60 51 ! to use 'getin' 61 !! **WRF a new stuff to be checked62 52 USE ioipsl_getincom 63 53 IMPLICIT NONE … … 71 61 #include "callkeys.h" 72 62 #include "surfdat.h" 73 #include "dimradmars.h" !!new74 #include "yomaer.h" !!new -- prob. pour remplir les proprietes63 #include "dimradmars.h" 64 #include "yomaer.h" 75 65 #include "datafile.h" 66 #ifdef MESOSCALE 76 67 #include "meso_slope.h" 77 #include "comsoil.h" !!**WRF -- needed to fill volcapa 78 68 #include "comsoil.h" !!MESOSCALE -- needed to fill volcapa 69 #include "meso_inc/meso_inc_inifisvar.F" 70 #endif 79 71 REAL prad,pg,pr,pcpp,pdaysec 80 81 INTEGER ngrid,nlayer,nq 82 83 REAL womeg,wmugaz,wdaysec 84 REAL wyear_day,wperiheli,waphelie,wperi_day,wobliquit 85 REAL wz0,wemin_turb,wlmixmin 86 REAL wemissiv,wemissiceN,wemissiceS,walbediceN,walbediceS 87 88 REAL wiceradiusN,wiceradiusS,wdtemisiceN,wdtemisiceS 89 REAL walbedodat(ngrid),wphisfi(ngrid) 90 REAL wzmea(ngrid),wzstd(ngrid),wzsig(ngrid) 91 REAL wzgam(ngrid),wzthe(ngrid) 92 REAL wtheta(ngrid),wpsi(ngrid) 93 REAL wvolcapa 94 REAL wdt 95 72 #ifndef MESOSCALE 73 REAL ptimestep 74 INTEGER day_ini 75 #endif 76 INTEGER ngrid,nlayer 96 77 REAL plat(ngrid),plon(ngrid),parea(ngridmx) 97 integer wday_ini98 78 INTEGER ig,ierr 99 100 INTEGER wecri_phys, wappel_phys101 79 102 80 ! EXTERNAL iniorbit,orbite … … 107 85 CHARACTER ch80*80 108 86 109 #ifdef MESOSCALE110 111 87 ! logical chem, h2o 112 88 … … 114 90 ! h2o = .false. 115 91 116 c ****WRF 117 c 118 c------------------------------------------------------ 119 c Fill some parameters in the 'include' files 120 c >> Do part of the job previously done by phyetat0.F 121 c >> Complete list of parameters is found in tabfi.F 122 c------------------------------------------------------ 123 c 124 c Values are defined in the module_model_constants.F WRF routine 125 c 126 ! in 'comcstfi.h' 127 rad=prad 92 rad=prad 128 93 cpp=pcpp 129 94 g=pg 130 95 r=pr 131 96 rcp=r/cpp 132 daysec=wdaysec 133 omeg=womeg 134 mugaz=wmugaz 135 print*,"check: rad,cpp,g,r,rcp,daysec,omeg,mugaz" 136 print*,rad,cpp,g,r,rcp,daysec,omeg,mugaz 137 138 ! in 'planet.h' 139 year_day=wyear_day 140 periheli=wperiheli 141 aphelie=waphelie 142 peri_day=wperi_day 143 obliquit=wobliquit 144 z0=wz0 145 emin_turb=wemin_turb 146 lmixmin=wlmixmin 147 print*,"check: year_day,periheli,aphelie,peri_day,obliquit" 148 print*,year_day,periheli,aphelie,peri_day,obliquit 149 print*,"check: z0,emin_turb,lmixmin" 150 print*,z0,emin_turb,lmixmin 151 152 ! in 'surfdat.h' 153 emissiv=wemissiv 154 emisice(1)=wemissiceN 155 emisice(2)=wemissiceS 156 albedice(1)=walbediceN 157 albedice(2)=walbediceS 158 iceradius(1)=wiceradiusN 159 iceradius(2)=wiceradiusS 160 dtemisice(1)=wdtemisiceN 161 dtemisice(2)=wdtemisiceS 162 print*,"check: emissiv,emisice,albedice,iceradius,dtemisice" 163 print*,emissiv,emisice,albedice,iceradius,dtemisice 164 165 c 166 c Values are defined in the WPS processing 167 c 168 albedodat(:)=walbedodat(:) 169 !!!!! ***WRF inertiedat was moved, new physics !! 170 !inertiedat(:)=winertiedat(:) 171 phisfi(:)=wphisfi(:) 172 print*,"check: albedodat(1),phisfi(1)" 173 print*,albedodat(1),phisfi(1) 174 print*,"check: albedodat(end),phisfi(end)" 175 print*,albedodat(ngrid),phisfi(ngrid) 176 177 ! NB: usually, gravity wave scheme is useless in mesoscale modeling 178 ! NB: we however keep the option for coarse grid case ... 179 zmea(:)=wzmea(:) 180 zstd(:)=wzstd(:) 181 zsig(:)=wzsig(:) 182 zgam(:)=wzgam(:) 183 zthe(:)=wzthe(:) 184 print*,"check: gw param" 185 print*,zmea(1),zmea(ngrid) 186 print*,zstd(1),zstd(ngrid) 187 print*,zsig(1),zsig(ngrid) 188 print*,zgam(1),zgam(ngrid) 189 print*,zthe(1),zthe(ngrid) 190 191 ! 192 ! in meso_slope.h 193 ! 194 theta_sl(:)=wtheta(:) 195 psi_sl(:)=wpsi(:) 196 print*,"check: theta_sl(1),psi_sl(1)" 197 print*,theta_sl(1),psi_sl(1) 198 print*,"check: theta_sl(end),psi_sl(end)" 199 print*,theta_sl(ngrid),psi_sl(ngrid) 200 201 ! 202 ! in comsoil.h 203 ! 204 volcapa=wvolcapa 205 print*,"check: volcapa" 206 print*,volcapa 207 208 209 c ****WRF 210 211 212 ! rad=prad 213 ! daysec=pdaysec 214 ! dtphys=ptimestep 215 ! cpp=pcpp 216 ! g=pg 217 ! r=pr 218 ! rcp=r/cpp 97 #ifndef MESOSCALE 98 daysec=pdaysec 99 dtphys=ptimestep 100 #else 101 #include "meso_inc/meso_inc_inifisini.F" 102 #endif 219 103 220 104 ! -------------------------------------------------------- … … 254 138 255 139 write(*,*) "Directory where external input files are:" 256 !datafile="/u/forget/WWW/datagcm/datafile" 257 ! NB: default 'datafile' is set in datafile.h 258 call getin("datadir",datafile) 140 #ifndef MESOSCALE 141 datafile="/u/forget/WWW/datagcm/datafile" 142 #endif 143 call getin("datadir",datafile) ! default path 259 144 write(*,*) " datafile = ",trim(datafile) 260 145 … … 283 168 284 169 write(*,*) "Save statistics in file stats.nc ?" 285 !callstats=.true. ! default value 170 #ifdef MESOSCALE 286 171 callstats=.false. ! default value 172 #else 173 callstats=.true. ! default value 174 #endif 287 175 call getin("callstats",callstats) 288 176 write(*,*) " callstats = ",callstats … … 639 527 8001 FORMAT(t5,a12,i8) 640 528 641 642 643 c ****WRF644 c*****************************************************645 c Since it comes from WRF settings, we have to646 c fill dtphys in the include file647 c It must be set now, because it is used afterwards648 c649 c Opportunity is taken to fill ecri_phys as well650 c*****************************************************651 dtphys=wdt*float(wappel_phys)652 print*,'Physical timestep (s) ',dtphys653 ecri_phys=wecri_phys654 print*,'Physical frequency for writing ',ecri_phys655 c656 c ****WRF657 658 659 529 PRINT* 660 530 PRINT*,'inifis: daysec',daysec … … 734 604 ! ------------------------ 735 605 736 ! in 'comgeomfi.h' 606 ! in 'comgeomfi.h' 737 607 CALL SCOPY(ngrid,plon,1,long,1) 738 608 CALL SCOPY(ngrid,plat,1,lati,1) … … 824 694 ! 825 695 ! RETURN 826 827 #endif828 829 696 END -
trunk/LMDZ.MARS/libf/phymars/meso_param_slope.F90
r47 r226 88 88 deg2rad = pi/180. 89 89 if ((theta_s > 90.) .or. (theta_s < 0.)) then 90 91 90 print *, 'please set theta_s between 0 and 90', theta_s 91 stop 92 92 endif 93 93 … … 162 162 ! 163 163 if (csza .ge. 0.5) then 164 165 166 167 168 169 else 170 171 172 173 164 mat_M(:,1) = (/ -0.264, 1.309, 0.208, -0.828 /) 165 mat_M(:,2) = (/ 1.291*sigma_s, -1.371*sigma_s, -0.581, 1.641 /) 166 mat_N(:,1) = (/ 0.911, -0.777, -0.223, 0.623 /) 167 mat_N(:,2) = (/ -0.933*sigma_s, 0.822*sigma_s, 0.514, -1.195 /) 168 169 else 170 mat_M(:,1) = (/ -0.373, 0.792, -0.095, 0.398 /) 171 mat_M(:,2) = (/ 1.389*sigma_s, -0.794*sigma_s, -0.325, 0.183 /) 172 mat_N(:,1) = (/ 1.079, 0.275, 0.419, -1.855 /) 173 mat_N(:,2) = (/ -1.076*sigma_s, -0.357*sigma_s, -0.075, 1.844 /) 174 174 endif 175 175 ! … … 184 184 ! low angles 185 185 ! 186 187 186 s_vector = (/ 1., exp(-taudust) , sin(0.0872664626), sin(0.0872664626)*exp(-taudust) /) 187 ratio = DOT_PRODUCT ( MATMUL( s_vector, mat_T), g_vector ) 188 188 ratio = 1. + (ratio - 1.)*deg2rad*theta_s/0.0872664626 189 189 else … … 191 191 ! general case 192 192 ! 193 193 ratio= DOT_PRODUCT ( MATMUL( s_vector, mat_T), g_vector ) 194 194 ! 195 195 ! NB: ratio= DOT_PRODUCT ( s_vector, MATMUL( mat_T, g_vector ) ) is equivalent -
trunk/LMDZ.MARS/libf/phymars/meso_physiq.F
r185 r226 1 SUBROUTINE meso_physiq(ngrid,nlayer,nq, 2 $ firstcall,lastcall, 3 $ wday_ini, 4 $ pday,ptime,ptimestep, 5 $ pplev,pplay,pphi, 6 $ pu,pv,pt,pq,pw, 7 $ wtnom, 8 $ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn, 9 $ wtsurf,wtsoil,wemis,wq2,wqsurf,wco2ice, 10 $ wisoil,wdsoil, 11 $ wecritphys, 1 SUBROUTINE meso_physiq( 2 $ ngrid,nlayer,nq 3 $ ,firstcall,lastcall 4 $ ,pday,ptime,ptimestep 5 $ ,pplev,pplay,pphi 6 $ ,pu,pv,pt,pq 7 $ ,pw 8 $ ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn 12 9 #ifdef MESOSCALE 13 $ output_tab2d, output_tab3d, 14 #endif 15 $ flag_LES) 16 17 10 #include "meso_inc/meso_inc_invar.F" 11 #endif 12 $ ) 18 13 19 14 IMPLICIT NONE 20 c=======================================================================21 c22 c CAREFUL: THIS IS A VERSION TO BE USED WITH WRF !!!23 c24 c ... CHECK THE ****WRF lines25 c26 15 c======================================================================= 27 16 c … … 75 64 c Radiatively active tracers: J.-B. Madeleine (10/2008-06/2009) 76 65 c Nb: See callradite.F for more information. 77 c new WRF version: Aymeric Spiga (01/2009)66 c Mesoscale version: Aymeric Spiga (2007 - 2011) -- check MESOSCALE flags 78 67 c 79 c80 68 c arguments: 81 69 c ---------- … … 107 95 c pw(ngrid,?) vertical velocity 108 96 c 109 c110 c ****WRF111 c day_ini,tsurf,tsoil,emis,q2,qsurf,co2ice are inputs112 c and locally saved variables113 c (no need to call phyetat0)114 c115 c116 97 c output: 117 98 c ------- … … 134 115 #include "comgeomfi.h" 135 116 #include "surfdat.h" 136 #include "comsoil.h" !!! new soil common117 #include "comsoil.h" 137 118 #include "comdiurn.h" 138 119 #include "callkeys.h" … … 154 135 #include "netcdf.inc" 155 136 156 !!!!**** SPECIFIC TO MESOSCALE157 137 #ifdef MESOSCALE 158 138 #include "meso_slope.h" 159 139 #include "wrf_output_2d.h" 160 140 #include "wrf_output_3d.h" 161 #endif162 163 141 #include "advtrac.h" !!! this is necessary for tracers (in dyn3d) 142 #include "meso_inc/meso_inc_var.F" 143 #endif 164 144 165 145 c Arguments : … … 177 157 REAL zh(ngridmx,nlayermx) ! potential temperature (K) 178 158 LOGICAL firstcall,lastcall 179 !!! ****WRF WRF specific to mesoscale 180 INTEGER wday_ini 181 REAL wtsurf(ngridmx) ! input only ay firstcall - output 182 REAL wtsoil(ngridmx,nsoilmx) 183 REAL wisoil(ngridmx,nsoilmx) !! new soil scheme 184 REAL wdsoil(ngridmx,nsoilmx) !! new soil scheme 185 REAL wco2ice(ngridmx) 186 REAL wemis(ngridmx) 187 REAL wqsurf(ngridmx,nqmx) 188 REAL wq2(ngridmx,nlayermx+1) 189 REAL wecritphys 190 #ifdef MESOSCALE 191 REAL output_tab2d(ngridmx,n2d) 192 REAL output_tab3d(ngridmx,nlayer,n3d) 193 #endif 194 REAL sl_ls, sl_lct, sl_lat, sl_tau, sl_alb, sl_the, sl_psi 195 REAL sl_fl0, sl_flu 196 REAL sl_ra, sl_di0 197 REAL sky 198 REAL hfx(ngridmx) !! pour LES avec isfflx!=0 199 REAL ust(ngridmx) !! pour LES avec isfflx!=0 200 LOGICAL flag_LES !! pour LES avec isfflx!=0 201 REAL qsurfice(ngridmx) !! pour diagnostics 202 real alpha,lay1 ! coefficients for building layers 203 integer iloop 204 INTEGER tracerset !!! this corresponds to config%mars 205 !!! ****WRF WRF specific to mesoscale 159 206 160 REAL pday 207 161 REAL ptime 208 162 logical tracerdyn 209 CHARACTER (len=20) :: wtnom(nqmx) ! tracer name210 163 211 164 c outputs: … … 371 324 372 325 c======================================================================= 373 #ifdef MESOSCALE374 326 375 327 c 1. Initialisation: … … 388 340 c read startfi 389 341 c ~~~~~~~~~~~~ 390 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 391 c ****WRF 392 c 393 c No need to use startfi.nc 394 c > part of the job of phyetat0 is done in inifis 395 c > remaining initializations are passed here from the WRF variables 396 c > beware, some operations were done by phyetat0 (ex: tracers) 397 c > if any problems, look in phyetat0 398 c 399 tsurf(:)=wtsurf(:) 400 PRINT*,'check: tsurf ',tsurf(1),tsurf(ngridmx) 401 tsoil(:,:)=wtsoil(:,:) 402 PRINT*,'check: tsoil ',tsoil(1,1),tsoil(ngridmx,nsoilmx) 403 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 404 !!!new physics 405 inertiedat(:,:)=wisoil(:,:) 406 PRINT*,'check: inert ',inertiedat(1,1),inertiedat(ngridmx,nsoilmx) 407 mlayer(0:nsoilmx-1)=wdsoil(1,:) 408 PRINT*,'check: midlayer ', mlayer(:) 409 !!!!!!!!!!!!!!!!! DONE in soil_setting.F 410 ! 1.5 Build layer(); following the same law as mlayer() 411 ! Assuming layer distribution follows mid-layer law: 412 ! layer(k)=lay1*alpha**(k-1) 413 lay1=sqrt(mlayer(0)*mlayer(1)) 414 alpha=mlayer(1)/mlayer(0) 415 do iloop=1,nsoilmx 416 layer(iloop)=lay1*(alpha**(iloop-1)) 417 enddo 418 419 PRINT*,'check: layer ', layer(:) 420 421 !!!!!!!!!!!!!!!!! DONE in soil_setting.F 422 tnom(:)=wtnom(:) !! est rempli dans advtrac.h 423 PRINT*,'check: tracernames ', tnom 424 !!!new physics 425 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 426 emis(:)=wemis(:) 427 PRINT*,'check: emis ',emis(1),emis(ngridmx) 428 q2(:,:)=wq2(:,:) 429 PRINT*,'check: q2 ',q2(1,1),q2(ngridmx,nlayermx+1) 430 qsurf(:,:)=wqsurf(:,:) 431 PRINT*,'check: qsurf ',qsurf(1,1),qsurf(ngridmx,nqmx) 432 co2ice(:)=wco2ice(:) 433 PRINT*,'check: co2 ',co2ice(1),co2ice(ngridmx) 434 day_ini=wday_ini 435 436 c artificially filling dyn3d/control.h is also required 437 c > iphysiq is put in WRF to be set easily (cf ptimestep) 438 c > day_step is simply deduced: 439 c 440 day_step=daysec/ptimestep 441 PRINT*,'Call to LMD physics:',day_step,' per Martian day' 442 c 443 iphysiq=ptimestep 444 c 445 ecritphy=wecritphys 446 PRINT*,'Write LMD physics each:',ecritphy,' seconds' 447 !!PRINT*,ecri_phys 448 !!PRINT*,float(ecri_phys) ... 449 !!renvoient tous deux des nombres absurdes 450 !!pourtant callkeys.h est inclus ... 451 !! 452 !!donc ecritphys est passe en argument ... 453 PRINT*,'Write LMD physics each:',ecritphy,' seconds' 454 c 455 !DO iq=1, nq 456 ! PRINT*, tnom(iq), pq(:,:,iq) 457 !ENDDO 458 459 c 460 c ****WRF 461 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 462 463 464 465 !! Read netcdf initial physical parameters. 466 ! CALL phyetat0 ("startfi.nc",0,0, 467 ! & nsoilmx,nq, 468 ! & day_ini,time_phys, 469 ! & tsurf,tsoil,emis,q2,qsurf,co2ice) 342 #ifndef MESOSCALE 343 ! Read netcdf initial physical parameters. 344 CALL phyetat0 ("startfi.nc",0,0, 345 & nsoilmx,nq, 346 & day_ini,time_phys, 347 & tsurf,tsoil,emis,q2,qsurf,co2ice) 348 #else 349 #include "meso_inc/meso_inc_ini.F" 350 #endif 470 351 471 352 if (pday.ne.day_ini) then … … 507 388 ENDIF ! end tracer 508 389 509 !!!!!! WRF WRF WRF MARS MARS 510 !!!!!! TEST TEST TEST TEST AS+JBM 28/02/11 511 !!!!!! TEST TEST TEST TEST AS+JBM 28/02/11 512 !!!!!! TEST TEST TEST TEST AS+JBM 28/02/11 513 !!!! 514 !!!! principe: une option 'caps=T' specifique au mesoscale 515 !!!! ... en vue d'un meso_initracer ???? 516 !!!! 517 !!!! depots permanents => albedo TES du PDS 518 !!!! depots saisonniers => alb_surfice (~0.4, cf plus bas) 519 !!!! [!!!! y compris pour les depots saisonniers sur les depots permanents] 520 !!!! 521 !!!! --> todo: il faut garder les depots saisonniers qui viennent 522 !!!! du GCM lorsqu'ils sont consequents 523 !!!! 524 IF ( caps .and. (igcm_h2o_ice .ne. 0) ) THEN 525 PRINT *, 'OVERWRITING watercaptag DEFINITION in INITRACER' 526 PRINT *, 'lat>70 et alb>0.26 => watercaptag=T' 527 !! Perennial H20 north cap defined by watercaptag=true (allows surface to be 528 !! hollowed by sublimation in vdifc). 529 do ig=1,ngridmx 530 qsurf(ig,igcm_h2o_ice)=0. !! on jette les inputs GCM 531 if ( (lati(ig)*180./pi.gt.70.) .and. 532 . (albedodat(ig).ge.0.26) ) then 533 watercaptag(ig)=.true. 534 dryness(ig) = 1. 535 else 536 watercaptag(ig)=.false. 537 dryness(ig) = 1. 538 endif ! (lati, albedodat) 539 end do ! (ngridmx) 540 ELSE ! (caps) 541 print *,'Blork !!!' 542 print *,'caps=T avec water=F ????' 543 ENDIF ! (caps) 544 !!!!!! TEST TEST TEST TEST AS+JBM 28/02/11 545 !!!!!! TEST TEST TEST TEST AS+JBM 28/02/11 546 !!!!!! TEST TEST TEST TEST AS+JBM 28/02/11 547 548 549 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 550 c ****WRF 551 c 552 c nosense in mesoscale modeling 553 c 554 cc Determining gridpoint near Viking Lander 1 (used for diagnostic only) 555 cc ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 556 c 557 c if(ngrid.ne.1) then 558 c latvl1= 22.27 559 c lonvl1= -47.94 560 c ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.) -2 )*iim + 561 c & int(1.5+(lonvl1+180)*iim/360.) 562 c write(*,*) 'Viking Lander 1 GCM point: lat,lon', 563 c & lati(ig_vl1)*180/pi, long(ig_vl1)*180/pi 564 c end if 565 c ****WRF 566 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 567 568 !!! 569 !!! WRF WRF WRF commented for smaller executables 570 !!! 571 !c Initialize thermospheric parameters 572 !c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 573 ! 574 ! if (callthermos) call param_read 575 576 390 #ifdef MESOSCALE 391 #include "meso_inc/meso_inc_caps.F" 392 #endif 393 394 #ifndef MESOSCALE 395 c Determining gridpoint near Viking Lander 1 (used for diagnostic only) 396 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 397 398 if(ngrid.ne.1) then 399 latvl1= 22.27 400 lonvl1= -47.94 401 ig_vl1= 1+ int( (1.5-(latvl1-90.)*jjm/180.) -2 )*iim + 402 & int(1.5+(lonvl1+180)*iim/360.) 403 write(*,*) 'Viking Lander 1 GCM point: lat,lon', 404 & lati(ig_vl1)*180/pi, long(ig_vl1)*180/pi 405 end if 406 #endif 407 408 #ifndef MESOSCALE 409 c Initialize thermospheric parameters 410 c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 411 412 if (callthermos) call param_read 413 #endif 577 414 c Initialize R and Cp as constant 578 415 … … 661 498 ENDDO 662 499 663 !!! 664 !!! WRF WRF WRF commented for smaller executables 665 !!! 666 !c----------------------------------------------------------------------- 667 !c 1.2.5 Compute mean mass, cp, and R 668 !c -------------------------------- 669 ! 670 ! if(photochem.or.callthermos) then 671 ! call concentrations(pplay,pt,pdt,pq,pdq,ptimestep) 672 ! endif 673 500 #ifndef MESOSCALE 501 c----------------------------------------------------------------------- 502 c 1.2.5 Compute mean mass, cp, and R 503 c -------------------------------- 504 505 if(photochem.or.callthermos) then 506 call concentrations(pplay,pt,pdt,pq,pdq,ptimestep) 507 endif 508 #endif 674 509 c----------------------------------------------------------------------- 675 510 c 2. Compute radiative tendencies : … … 680 515 IF( MOD(icount-1,iradia).EQ.0) THEN 681 516 517 #ifdef MESOSCALE 682 518 write (*,*) 'call radiative transfer' 683 519 #endif 684 520 c Local Solar zenith angle 685 521 c ~~~~~~~~~~~~~~~~~~~~~~~~ … … 716 552 c Radiative transfer 717 553 c ------------------ 718 cc 719 cc **WRF: desormais dust_opacity est dans callradite -- modifications 720 cc nveaux arguments: tauref,tau,aerosol,rice,nuice 721 cc 554 722 555 CALL callradite(icount,ngrid,nlayer,nq,zday,zls,pq,albedo, 723 556 $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, … … 725 558 & tauref,tau,aerosol,ccn,rdust,rice,nuice) 726 559 727 c write(*,*) icount,ngrid,nlayer,nq,zday,zls,pq,albedo, 728 c $ emis,mu0,pplev,pplay,pt,tsurf,fract,dist_sol,igout, 729 c $ zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,fluxtop_sw, 730 c & tauref,tau,aerosol,rice,nuice 731 c write(*,*) fluxsurf_lw 732 733 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 734 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 735 ccccc 736 ccccc PARAM SLOPE : Insolation (direct + scattered) 737 ccccc 738 DO ig=1,ngrid 739 sl_the = theta_sl(ig) 740 IF (sl_the .ne. 0.) THEN 741 ztim1=fluxsurf_sw(ig,1)+fluxsurf_sw(ig,2) 742 DO l=1,2 743 sl_lct = ptime*24. + 180.*long(ig)/pi/15. 744 sl_ra = pi*(1.0-sl_lct/12.) 745 sl_lat = 180.*lati(ig)/pi 746 sl_tau = tau(ig,1) 747 sl_alb = albedo(ig,l) 748 sl_psi = psi_sl(ig) 749 sl_fl0 = fluxsurf_sw(ig,l) 750 sl_di0 = 0. 751 if (mu0(ig) .gt. 0.) then 752 sl_di0 = mu0(ig)*(exp(-sl_tau/mu0(ig))) 753 sl_di0 = sl_di0*1370./dist_sol/dist_sol 754 sl_di0 = sl_di0/ztim1 755 sl_di0 = fluxsurf_sw(ig,l)*sl_di0 756 endif 757 ! sait-on jamais (a cause des arrondis) 758 if (sl_fl0 .lt. sl_di0) sl_di0=sl_fl0 759 !!!!!!!!!!!!!!!!!!!!!!!!!! 760 CALL meso_param_slope( mu0(ig), declin, sl_ra, sl_lat, 761 & sl_tau, sl_alb, 762 & sl_the, sl_psi, sl_di0, sl_fl0, sl_flu) 763 !!!!!!!!!!!!!!!!!!!!!!!!!! 764 fluxsurf_sw(ig,l) = sl_flu 765 !! sl_ls = 180.*zls/pi 766 !! sl_lct = ptime*24. + 180.*long(ig)/pi/15. 767 !! sl_lat = 180.*lati(ig)/pi 768 !! sl_tau = tau(ig,1) 769 !! sl_alb = albedo(ig,l) 770 !! sl_the = theta_sl(ig) 771 !! sl_psi = psi_sl(ig) 772 !! sl_fl0 = fluxsurf_sw(ig,l) 773 !! CALL param_slope_full(sl_ls, sl_lct, sl_lat, 774 !! & sl_tau, sl_alb, 775 !! & sl_the, sl_psi, sl_fl0, sl_flu) 776 ENDDO 777 !!! compute correction on IR flux as well 778 sky= (1.+cos(pi*theta_sl(ig)/180.))/2. 779 fluxsurf_lw(ig)= fluxsurf_lw(ig)*sky 780 ENDIF 781 ENDDO 782 ccccc 783 ccccc PARAM SLOPE 784 ccccc 785 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 786 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 787 560 #ifdef MESOSCALE 561 #include "meso_inc/meso_inc_slope.F" 562 #endif 788 563 789 564 c CO2 near infrared absorption … … 801 576 $ +fluxsurf_sw(ig,1)*(1.-albedo(ig,1)) 802 577 $ +fluxsurf_sw(ig,2)*(1.-albedo(ig,2)) 803 804 !print*,'RAD ', fluxrad_sky(ig)805 !print*,'LW ', emis(ig)*fluxsurf_lw(ig)806 !print*,'SW ', fluxsurf_sw(ig,2)*(1.-albedo(ig,2))807 808 578 ENDDO 809 579 … … 823 593 ENDIF 824 594 825 826 827 595 ENDIF ! of if(mod(icount-1,iradia).eq.0) 828 596 … … 838 606 $ stephan*zplanck(ig)*zplanck(ig) 839 607 fluxrad(ig)=fluxrad_sky(ig)-zplanck(ig) 840 cccc 841 cccc param slope 842 cccc 608 #ifdef MESOSCALE 609 !!!! param slope 843 610 sky= (1.+cos(pi*theta_sl(ig)/180.))/2. 844 611 fluxrad(ig)=fluxrad(ig)+(1.-sky)*zplanck(ig) 845 cccc 846 cccc 847 cccc 612 #endif 848 613 ENDDO 849 850 614 851 615 DO l=1,nlayer … … 857 621 ENDIF ! of IF (callrad) 858 622 859 !!! 860 !!! WRF WRF WRF commented for smaller executables 861 !!! 862 !c----------------------------------------------------------------------- 863 !c 3. Gravity wave and subgrid scale topography drag : 864 !c ------------------------------------------------- 865 ! 866 ! 867 ! IF(calllott)THEN 868 ! 869 ! CALL calldrag_noro(ngrid,nlayer,ptimestep, 870 ! & pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw) 871 ! 872 ! DO l=1,nlayer 873 ! DO ig=1,ngrid 874 ! pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l) 875 ! pdu(ig,l)=pdu(ig,l)+zdugw(ig,l) 876 ! pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l) 877 ! ENDDO 878 ! ENDDO 879 ! ENDIF 880 623 #ifndef MESOSCALE 624 c----------------------------------------------------------------------- 625 c 3. Gravity wave and subgrid scale topography drag : 626 c ------------------------------------------------- 627 628 629 IF(calllott)THEN 630 631 CALL calldrag_noro(ngrid,nlayer,ptimestep, 632 & pplay,pplev,pt,pu,pv,zdtgw,zdugw,zdvgw) 633 634 DO l=1,nlayer 635 DO ig=1,ngrid 636 pdv(ig,l)=pdv(ig,l)+zdvgw(ig,l) 637 pdu(ig,l)=pdu(ig,l)+zdugw(ig,l) 638 pdt(ig,l)=pdt(ig,l)+zdtgw(ig,l) 639 ENDDO 640 ENDDO 641 ENDIF 642 #endif 881 643 c----------------------------------------------------------------------- 882 644 c 4. Vertical diffusion (turbulent mixing): 883 645 c ----------------------------------------- 884 c 646 885 647 IF (calldifv) THEN 886 887 648 888 649 DO ig=1,ngrid 889 650 zflubid(ig)=fluxrad(ig)+fluxgrd(ig) 890 !write (*,*), fluxrad(ig), fluxgrd(ig), zflubid(ig)891 651 ENDDO 892 652 … … 898 658 enddo 899 659 enddo 900 660 901 661 c Calling vdif (Martian version WITH CO2 condensation) 902 662 CALL vdifc(ngrid,nlayer,nq,co2ice,zpopsk, … … 908 668 & zdqdif,zdqsdif) 909 669 910 DO ig=1,ngrid 911 !! sensible heat flux in W/m2 912 hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig) 913 !! u star in similarity theory in m/s 914 ust(ig) = 0.4 915 . * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ) 916 . / log( 1.E+0 + zzlay(ig,1)/z0 ) 917 ENDDO 918 919 ! write (*,*) 'PHYS HFX cp zdts', hfx(100), zflubid(100), 920 ! . capcal(100), 921 ! . zdtsdif(100) 922 ! write (*,*) 'PHYS UST', ust(100) 923 924 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 925 !!! LES LES 926 IF (flag_LES) THEN 927 928 write (*,*) '************************************************' 929 write (*,*) '** LES mode: the difv part is only used to' 930 write (*,*) '** provide HFX and UST to the dynamics' 931 write (*,*) '** NB: - dudif, dvdif, dhdif, dqdif are set to 0' 932 write (*,*) '** - tsurf is updated' 933 write (*,*) '************************************************' 934 935 !!! 936 !!! WRF WRF LES LES : en fait le subgrid scale n'etait pas mis a zero !! 937 !!! 938 DO ig=1,ngrid 939 ! !! sensible heat flux in W/m2 940 ! hfx(ig) = zflubid(ig)-capcal(ig)*zdtsdif(ig) 941 ! !! u star in similarity theory in m/s 942 ! ust(ig) = 0.4 943 ! . * sqrt( pu(ig,1)*pu(ig,1) + pv(ig,1)*pv(ig,1) ) 944 ! . / log( 1.E+0 + zzlay(ig,1)/z0 ) 945 ! 946 DO l=1,nlayer 947 zdvdif(ig,l) = 0. 948 zdudif(ig,l) = 0. 949 zdhdif(ig,l) = 0. 950 DO iq=1, nq 951 zdqdif(ig,l,iq) = 0. 952 zdqsdif(ig,iq) = 0. !! sortir de la boucle 953 ENDDO 954 ENDDO 955 ! 956 ENDDO 957 !write (*,*) 'RAD ',fluxrad(igout) 958 !write (*,*) 'GRD ',fluxgrd(igout) 959 !write (*,*) 'dTs/dt ',capcal(igout)*zdtsurf(igout) 960 !write (*,*) 'HFX ', hfx(igout) 961 !write (*,*) 'UST ', ust(igout) 962 ENDIF 963 !!! LES LES 964 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 965 670 #ifdef MESOSCALE 671 #include "meso_inc/meso_inc_les.F" 672 #endif 966 673 DO l=1,nlayer 967 674 DO ig=1,ngrid … … 975 682 ENDDO 976 683 977 DO ig=1,ngrid978 zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig)979 ENDDO684 DO ig=1,ngrid 685 zdtsurf(ig)=zdtsurf(ig)+zdtsdif(ig) 686 ENDDO 980 687 981 688 if (tracer) then … … 999 706 s (fluxrad(ig)+fluxgrd(ig))/capcal(ig) 1000 707 ENDDO 1001 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 708 #ifdef MESOSCALE 1002 709 IF (flag_LES) THEN 1003 710 write(*,*) 'LES mode !' … … 1005 712 STOP 1006 713 ENDIF 1007 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 714 #endif 1008 715 ENDIF ! of IF (calldifv) 1009 716 … … 1020 727 $ pplay,pplev,pphi,zpopsk, 1021 728 $ pdu_th,pdv_th,pdt_th,pdq_th,lmax_th,zmax_th, 1022 $ dtke_th,hfmax_th,wmax_th) 1023 729 $ dtke_th,hfmax_th,wmax_th) 730 1024 731 DO l=1,nlayer 1025 732 DO ig=1,ngrid … … 1103 810 $ co2ice,albedo,emis, 1104 811 $ zdtc,zdtsurfc,pdpsrf,zduc,zdvc,zdqc, 1105 $ 812 $ fluxsurf_sw,zls) 1106 813 1107 814 DO l=1,nlayer … … 1178 885 c ------------------ 1179 886 1180 !!! 1181 !!! WRF WRF WRF commented for smaller executables 1182 !!! 1183 !c -------------- 1184 !c photochemistry : 1185 !c -------------- 1186 ! IF (photochem .or. thermochem) then 1187 ! call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, 1188 ! & zzlay,zday,pq,pdq,rice, 1189 ! & zdqchim,zdqschim,zdqcloud,zdqscloud) 1190 !!NB: Photochemistry includes condensation of H2O2 1191 ! 1192 ! ! increment values of tracers: 1193 ! DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry 1194 ! ! tracers is zero anyways 1195 ! DO l=1,nlayer 1196 ! DO ig=1,ngrid 1197 ! pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq) 1198 ! ENDDO 1199 ! ENDDO 1200 ! ENDDO ! of DO iq=1,nq 1201 ! ! add condensation tendency for H2O2 1202 ! if (igcm_h2o2.ne.0) then 1203 ! DO l=1,nlayer 1204 ! DO ig=1,ngrid 1205 ! pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2) 1206 ! & +zdqcloud(ig,l,igcm_h2o2) 1207 ! ENDDO 1208 ! ENDDO 1209 ! endif 1210 ! 1211 ! ! increment surface values of tracers: 1212 ! DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry 1213 ! ! tracers is zero anyways 1214 ! DO ig=1,ngrid 1215 ! dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq) 1216 ! ENDDO 1217 ! ENDDO ! of DO iq=1,nq 1218 ! ! add condensation tendency for H2O2 1219 ! if (igcm_h2o2.ne.0) then 1220 ! DO ig=1,ngrid 1221 ! dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2) 1222 ! & +zdqscloud(ig,igcm_h2o2) 1223 ! ENDDO 1224 ! endif 1225 ! 1226 ! END IF ! of IF (photochem.or.thermochem) 887 #ifndef MESOSCALE 888 c -------------- 889 c photochemistry : 890 c -------------- 891 IF (photochem .or. thermochem) then 892 call calchim(ptimestep,pplay,pplev,pt,pdt,dist_sol,mu0, 893 & zzlay,zday,pq,pdq,rice, 894 & zdqchim,zdqschim,zdqcloud,zdqscloud) 895 !NB: Photochemistry includes condensation of H2O2 896 897 ! increment values of tracers: 898 DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry 899 ! tracers is zero anyways 900 DO l=1,nlayer 901 DO ig=1,ngrid 902 pdq(ig,l,iq)=pdq(ig,l,iq)+zdqchim(ig,l,iq) 903 ENDDO 904 ENDDO 905 ENDDO ! of DO iq=1,nq 906 ! add condensation tendency for H2O2 907 if (igcm_h2o2.ne.0) then 908 DO l=1,nlayer 909 DO ig=1,ngrid 910 pdq(ig,l,igcm_h2o2)=pdq(ig,l,igcm_h2o2) 911 & +zdqcloud(ig,l,igcm_h2o2) 912 ENDDO 913 ENDDO 914 endif 915 916 ! increment surface values of tracers: 917 DO iq=1,nq ! loop on all tracers; tendencies for non-chemistry 918 ! tracers is zero anyways 919 DO ig=1,ngrid 920 dqsurf(ig,iq)=dqsurf(ig,iq)+zdqschim(ig,iq) 921 ENDDO 922 ENDDO ! of DO iq=1,nq 923 ! add condensation tendency for H2O2 924 if (igcm_h2o2.ne.0) then 925 DO ig=1,ngrid 926 dqsurf(ig,igcm_h2o2)=dqsurf(ig,igcm_h2o2) 927 & +zdqscloud(ig,igcm_h2o2) 928 ENDDO 929 endif 930 931 END IF ! of IF (photochem.or.thermochem) 932 #endif 1227 933 1228 934 c 7c. Aerosol particles … … 1296 1002 endif ! of if (tracer) 1297 1003 1298 !!! 1299 !!! WRF WRF WRF commented for smaller executables 1300 !!! 1301 !c----------------------------------------------------------------------- 1302 !c 8. THERMOSPHERE CALCULATION 1303 !c----------------------------------------------------------------------- 1304 ! 1305 ! if (callthermos) then 1306 ! call thermosphere(pplev,pplay,dist_sol, 1307 ! $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay, 1308 ! & pt,pq,pu,pv,pdt,pdq, 1309 ! $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff) 1310 ! 1311 ! DO l=1,nlayer 1312 ! DO ig=1,ngrid 1313 ! dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l) 1314 ! pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l) 1315 ! & +zdteuv(ig,l) 1316 ! pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l) 1317 ! pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l) 1318 ! DO iq=1, nq 1319 ! pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq) 1320 ! ENDDO 1321 ! ENDDO 1322 ! ENDDO 1323 ! 1324 ! endif ! of if (callthermos) 1004 #ifndef MESOSCALE 1005 c----------------------------------------------------------------------- 1006 c 8. THERMOSPHERE CALCULATION 1007 c----------------------------------------------------------------------- 1008 1009 if (callthermos) then 1010 call thermosphere(pplev,pplay,dist_sol, 1011 $ mu0,ptimestep,ptime,zday,tsurf,zzlev,zzlay, 1012 & pt,pq,pu,pv,pdt,pdq, 1013 $ zdteuv,zdtconduc,zdumolvis,zdvmolvis,zdqmoldiff) 1014 1015 DO l=1,nlayer 1016 DO ig=1,ngrid 1017 dtrad(ig,l)=dtrad(ig,l)+zdteuv(ig,l) 1018 pdt(ig,l)=pdt(ig,l)+zdtconduc(ig,l) 1019 & +zdteuv(ig,l) 1020 pdv(ig,l)=pdv(ig,l)+zdvmolvis(ig,l) 1021 pdu(ig,l)=pdu(ig,l)+zdumolvis(ig,l) 1022 DO iq=1, nq 1023 pdq(ig,l,iq)=pdq(ig,l,iq)+zdqmoldiff(ig,l,iq) 1024 ENDDO 1025 ENDDO 1026 ENDDO 1027 1028 endif ! of if (callthermos) 1029 #endif 1325 1030 1326 1031 c----------------------------------------------------------------------- … … 1341 1046 1342 1047 IF (tracer.AND.water.AND.(ngridmx.NE.1)) THEN 1343 ! if (caps.and.(obliquit.lt.27.)) then 1344 ! ! NB: Updated surface pressure, at grid point 'ngrid', is 1345 ! ! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep 1346 ! tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095* 1347 ! & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) 1348 ! endif 1349 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1350 !!!!! note WRF MESOSCALE AYMERIC -- mot cle "caps" 1351 !!!!! watercaptag n'est plus utilise que dans vdifc 1352 !!!!! ... pour que la sublimation ne soit pas stoppee 1353 !!!!! ... dans la calotte permanente nord si qsurf=0 1354 !!!!! on desire garder cet effet regle par caps=T 1355 !!!!! on a donc commente "if (caps.and.(obliquit.lt.27.))" ci-dessus 1356 !!!!! --- remplacer ces lignes par qqch de plus approprie 1357 !!!!! si on s attaque a la calotte polaire sud 1358 !!!!! pas d'autre occurrence majeure du mot-cle "caps" 1359 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 1360 1048 #ifndef MESOSCALE 1049 if (caps.and.(obliquit.lt.27.)) then 1050 ! NB: Updated surface pressure, at grid point 'ngrid', is 1051 ! ps(ngrid)=pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep 1052 tsurf(ngrid)=1./(1./136.27-r/5.9e+5*alog(0.0095* 1053 & (pplev(ngrid,1)+pdpsrf(ngrid)*ptimestep))) 1054 endif 1055 #endif 1361 1056 c ------------------------------------------------------------- 1362 1057 c Change of surface albedo (set to 0.4) in case of ground frost … … 1365 1060 c ALWAYS PLACE these lines after newcondens !!! 1366 1061 c ------------------------------------------------------------- 1367 c **WRF : OK avec le mesoscale, pas d'indices bizarres au pole1368 1062 do ig=1,ngrid 1369 1063 if ((co2ice(ig).eq.0).and. … … 1429 1123 ENDDO 1430 1124 1431 c Compute surface stress : (NB: z0 is a common in planete.h)1125 c Compute surface stress : (NB: z0 is a common in surfdat.h) 1432 1126 c DO ig=1,ngrid 1433 c cd = (0.4/log(zzlay(ig,1)/z0 ))**21127 c cd = (0.4/log(zzlay(ig,1)/z0(ig)))**2 1434 1128 c zstress(ig) = rho(ig,1)*cd*(zu(ig,1)**2 + zv(ig,1)**2) 1435 1129 c ENDDO … … 1437 1131 c Sum of fluxes in solar spectral bands (for output only) 1438 1132 DO ig=1,ngrid 1439 1440 1133 fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2) 1134 fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2) 1441 1135 ENDDO 1442 1136 c ******* TEST ****************************************************** … … 1478 1172 1479 1173 IF (ngrid.NE.1) THEN 1480 print*,'Ls =',zls*180./pi, 1481 & ' tauref(700 Pa,lat=0) =',tauref(ngrid/2)!, 1482 ! & ' tau(Viking1) =',tau(ig_vl1,1) 1483 1484 1174 print*,'Ls =',zls*180./pi 1175 & ,' tauref(700 Pa,lat=0) =',tauref(ngrid/2) 1176 #ifndef MESOSCALE 1177 & ,' tau(Viking1) =',tau(ig_vl1,1) 1178 #endif 1179 1180 #ifndef MESOSCALE 1485 1181 c ------------------------------------------------------------------- 1486 1182 c Writing NetCDF file "RESTARTFI" at the end of the run … … 1492 1188 c thus we store for time=time+dtvr 1493 1189 1494 1495 !!! 1496 !!! WRF WRF WRF WRF 1497 !!! 1498 ! IF(lastcall) THEN 1499 ! ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) 1500 ! write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin 1501 ! call physdem1("restartfi.nc",long,lati,nsoilmx,nq, 1502 ! . ptimestep,pday, 1503 ! . ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf, 1504 ! . area,albedodat,inertiedat,zmea,zstd,zsig, 1505 ! . zgam,zthe) 1506 ! ENDIF 1507 1508 1190 IF(lastcall) THEN 1191 ztime_fin = ptime + ptimestep/(float(iphysiq)*daysec) 1192 write(*,*)'PHYSIQ: for physdem ztime_fin =',ztime_fin 1193 call physdem1("restartfi.nc",long,lati,nsoilmx,nq, 1194 . ptimestep,pday, 1195 . ztime_fin,tsurf,tsoil,co2ice,emis,q2,qsurf, 1196 . area,albedodat,inertiedat,zmea,zstd,zsig, 1197 . zgam,zthe) 1198 ENDIF 1199 #endif 1509 1200 1510 1201 c ------------------------------------------------------------------- … … 1557 1248 c "stats") only possible in 3D runs ! 1558 1249 1250 #ifndef MESOSCALE 1559 1251 IF (callstats) THEN 1560 1252 1561 write(*,*) 'callstats' 1562 1563 ! call wstats(ngrid,"ps","Surface pressure","Pa",2,ps) 1564 ! call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf) 1565 ! call wstats(ngrid,"co2ice","CO2 ice cover", 1566 ! & "kg.m-2",2,co2ice) 1567 ! call wstats(ngrid,"fluxsurf_lw", 1568 ! & "Thermal IR radiative flux to surface","W.m-2",2, 1569 ! & fluxsurf_lw) 1570 ! call wstats(ngrid,"fluxsurf_sw", 1571 ! & "Solar radiative flux to surface","W.m-2",2, 1572 ! & fluxsurf_sw_tot) 1573 ! call wstats(ngrid,"fluxtop_lw", 1574 ! & "Thermal IR radiative flux to space","W.m-2",2, 1575 ! & fluxtop_lw) 1576 ! call wstats(ngrid,"fluxtop_sw", 1577 ! & "Solar radiative flux to space","W.m-2",2, 1578 ! & fluxtop_sw_tot) 1579 ! call wstats(ngrid,"taudustvis", 1580 ! & "Dust optical depth"," ",2,tau(1,1)) 1581 ! call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) 1582 ! call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) 1583 ! call wstats(ngrid,"v","Meridional (North-South) wind", 1584 ! & "m.s-1",3,zv) 1585 !c call wstats(ngrid,"w","Vertical (down-up) wind", 1586 !c & "m.s-1",3,pw) 1587 ! call wstats(ngrid,"rho","Atmospheric density","none",3,rho) 1588 !c call wstats(ngrid,"pressure","Pressure","Pa",3,pplay) 1589 !c call wstats(ngrid,"q2", 1590 !c & "Boundary layer eddy kinetic energy", 1591 !c & "m2.s-2",3,q2) 1592 !c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2, 1593 !c & emis) 1594 !c call wstats(ngrid,"ssurf","Surface stress","N.m-2", 1595 !c & 2,zstress) 1596 ! 1597 ! if (tracer) then 1598 ! if (water) then 1599 ! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap) 1600 ! & *mugaz/mmol(igcm_h2o_vap) 1601 ! call wstats(ngrid,"vmr_h2ovapor", 1602 ! & "H2O vapor volume mixing ratio","mol/mol", 1603 ! & 3,vmr) 1604 ! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) 1605 ! & *mugaz/mmol(igcm_h2o_ice) 1606 ! call wstats(ngrid,"vmr_h2oice", 1607 ! & "H2O ice volume mixing ratio","mol/mol", 1608 ! & 3,vmr) 1609 ! 1610 ! call wstats(ngrid,"mtot", 1611 ! & "total mass of water vapor","kg/m2", 1612 ! & 2,mtot) 1613 ! call wstats(ngrid,"icetot", 1614 ! & "total mass of water ice","kg/m2", 1615 ! & 2,icetot) 1616 !c If activice is true, tauTES is computed in aeropacity.F; 1617 ! if (.not.activice) then 1618 ! call wstats(ngrid,"tauTES", 1619 ! & "tau abs 825 cm-1","", 1620 ! & 2,tauTES) 1621 ! endif ! of if (activice) 1622 ! 1623 ! endif ! of if (water) 1624 ! 1625 ! if (thermochem.or.photochem) then 1626 ! do iq=1,nq 1627 ! if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or. 1628 ! . (noms(iq).eq."co").or.(noms(iq).eq."n2").or. 1629 ! . (noms(iq).eq."h2").or. 1630 ! . (noms(iq).eq."o3")) then 1631 ! do l=1,nlayer 1632 ! do ig=1,ngrid 1633 ! vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq) 1634 ! end do 1635 ! end do 1636 ! call wstats(ngrid,"vmr_"//trim(noms(iq)), 1637 ! . "Volume mixing ratio","mol/mol",3,vmr) 1638 ! endif 1639 ! enddo 1640 ! endif ! of if (thermochem.or.photochem) 1641 ! 1642 ! endif ! of if (tracer) 1643 ! 1644 ! IF(lastcall) THEN 1645 ! write (*,*) "Writing stats..." 1646 ! call mkstats(ierr) 1647 ! ENDIF 1253 call wstats(ngrid,"ps","Surface pressure","Pa",2,ps) 1254 call wstats(ngrid,"tsurf","Surface temperature","K",2,tsurf) 1255 call wstats(ngrid,"co2ice","CO2 ice cover", 1256 & "kg.m-2",2,co2ice) 1257 call wstats(ngrid,"fluxsurf_lw", 1258 & "Thermal IR radiative flux to surface","W.m-2",2, 1259 & fluxsurf_lw) 1260 call wstats(ngrid,"fluxsurf_sw", 1261 & "Solar radiative flux to surface","W.m-2",2, 1262 & fluxsurf_sw_tot) 1263 call wstats(ngrid,"fluxtop_lw", 1264 & "Thermal IR radiative flux to space","W.m-2",2, 1265 & fluxtop_lw) 1266 call wstats(ngrid,"fluxtop_sw", 1267 & "Solar radiative flux to space","W.m-2",2, 1268 & fluxtop_sw_tot) 1269 call wstats(ngrid,"temp","Atmospheric temperature","K",3,zt) 1270 call wstats(ngrid,"u","Zonal (East-West) wind","m.s-1",3,zu) 1271 call wstats(ngrid,"v","Meridional (North-South) wind", 1272 & "m.s-1",3,zv) 1273 call wstats(ngrid,"w","Vertical (down-up) wind", 1274 & "m.s-1",3,pw) 1275 call wstats(ngrid,"rho","Atmospheric density","none",3,rho) 1276 call wstats(ngrid,"pressure","Pressure","Pa",3,pplay) 1277 c call wstats(ngrid,"q2", 1278 c & "Boundary layer eddy kinetic energy", 1279 c & "m2.s-2",3,q2) 1280 c call wstats(ngrid,"emis","Surface emissivity","w.m-1",2, 1281 c & emis) 1282 c call wstats(ngrid,"ssurf","Surface stress","N.m-2", 1283 c & 2,zstress) 1284 c call wstats(ngrid,"sw_htrt","sw heat.rate", 1285 c & "W.m-2",3,zdtsw) 1286 c call wstats(ngrid,"lw_htrt","lw heat.rate", 1287 c & "W.m-2",3,zdtlw) 1288 1289 if (tracer) then 1290 if (water) then 1291 vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_vap) 1292 & *mugaz/mmol(igcm_h2o_vap) 1293 call wstats(ngrid,"vmr_h2ovapor", 1294 & "H2O vapor volume mixing ratio","mol/mol", 1295 & 3,vmr) 1296 vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) 1297 & *mugaz/mmol(igcm_h2o_ice) 1298 call wstats(ngrid,"vmr_h2oice", 1299 & "H2O ice volume mixing ratio","mol/mol", 1300 & 3,vmr) 1301 call wstats(ngrid,"h2o_ice_s", 1302 & "surface h2o_ice","kg/m2", 1303 & 2,qsurf(1,igcm_h2o_ice)) 1304 1305 call wstats(ngrid,"mtot", 1306 & "total mass of water vapor","kg/m2", 1307 & 2,mtot) 1308 call wstats(ngrid,"icetot", 1309 & "total mass of water ice","kg/m2", 1310 & 2,icetot) 1311 call wstats(ngrid,"reffice", 1312 & "Mean reff","m", 1313 & 2,rave) 1314 c call wstats(ngrid,"rice", 1315 c & "Ice particle size","m", 1316 c & 3,rice) 1317 c If activice is true, tauTES is computed in aeropacity.F; 1318 if (.not.activice) then 1319 call wstats(ngrid,"tauTESap", 1320 & "tau abs 825 cm-1","", 1321 & 2,tauTES) 1322 endif 1323 1324 endif ! of if (water) 1325 1326 if (thermochem.or.photochem) then 1327 do iq=1,nq 1328 if ((noms(iq).eq."o").or.(noms(iq).eq."co2").or. 1329 . (noms(iq).eq."co").or.(noms(iq).eq."n2").or. 1330 . (noms(iq).eq."h2").or. 1331 . (noms(iq).eq."o3")) then 1332 do l=1,nlayer 1333 do ig=1,ngrid 1334 vmr(ig,l)=zq(ig,l,iq)*mmean(ig,l)/mmol(iq) 1335 end do 1336 end do 1337 call wstats(ngrid,"vmr_"//trim(noms(iq)), 1338 . "Volume mixing ratio","mol/mol",3,vmr) 1339 endif 1340 enddo 1341 endif ! of if (thermochem.or.photochem) 1342 1343 endif ! of if (tracer) 1344 1345 IF(lastcall) THEN 1346 write (*,*) "Writing stats..." 1347 call mkstats(ierr) 1348 ENDIF 1648 1349 1649 1350 ENDIF !if callstats 1351 #endif 1650 1352 1651 1353 c (Store EOF for Mars Climate database software) … … 1654 1356 ENDIF 1655 1357 1656 ccc**************** WRF OUTPUT ************************** 1657 ccc**************** WRF OUTPUT ************************** 1658 ccc**************** WRF OUTPUT ************************** 1659 !do ig=1,ngrid 1660 ! wtsurf(ig) = tsurf(ig) !! surface temperature 1661 ! wco2ice(ig) = co2ice(ig) !! co2 ice 1662 ! 1663 ! !!! specific to WRF WRF WRF 1664 ! !!! just to output water ice on surface 1665 ! !!! uncomment the Registry entry 1666 ! IF (igcm_h2o_ice .ne. 0) qsurfice(ig) = qsurf(ig,igcm_h2o_ice) 1667 ! 1668 ! !!! "VMR_ICE" "VOL. MIXING RATIO ICE" "ppm" 1669 ! IF (igcm_h2o_ice .ne. 0) THEN 1670 ! vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice)*mugaz/mmol(igcm_h2o_ice) 1671 ! ENDIF 1672 ! 1673 !enddo 1358 #ifdef MESOSCALE 1674 1359 wtsurf(1:ngrid) = tsurf(1:ngrid) !! surface temperature 1675 1360 wco2ice(1:ngrid) = co2ice(1:ngrid) !! co2 ice … … 1680 1365 . * mugaz / mmol(igcm_h2o_ice) 1681 1366 ENDIF 1682 1683 c 1684 c THIS INCLUDE IS AUTOMATICALLY GENERATED FROM REGISTRY 1685 c 1367 c AUTOMATICALLY GENERATED FROM REGISTRY 1686 1368 #include "fill_save.inc" 1687 c 1688 ccc**************** WRF OUTPUT ************************** 1689 ccc**************** WRF OUTPUT ************************** 1690 ccc**************** WRF OUTPUT ************************** 1691 1692 1693 cc----------------------------------- 1694 cc you can still use meso_WRITEDIAGFI (e.g. for debugging purpose), 1695 cc though this is not the default strategy now 1696 cc----------------------------------- 1697 cc please use cudt in namelist.input to set frequency of outputs 1698 cc----------------------------------- 1699 cc BEWARE: if at least one call to meso_WRITEDIAGFI is performed, 1700 cc cudt cannot be 0 - otherwise you'll get a "Floating exception" 1701 cc----------------------------------- 1702 ! call meso_WRITEDIAGFI(ngrid,"tauref", 1703 ! . "tauref","W.m-2",2, 1704 ! . tauref) 1705 ! call meso_WRITEDIAGFI(ngrid,"dtrad", 1706 ! . "dtrad","W.m-2",2, 1707 ! . dtrad) 1708 c call meso_WRITEDIAGFI(ngrid,"tsurf", 1709 c . "tsurf","K",2, 1710 c . tsurf) 1711 c 1712 ! call meso_WRITEDIAGFI(ngrid,"zt", 1713 ! . "zt","W.m-2",3, 1714 ! . zt) 1715 ! call meso_WRITEDIAGFI(ngrid,"zdtlw", 1716 ! . "zdtlw","W.m-2",2, 1717 ! . zdtlw) 1718 ! call meso_WRITEDIAGFI(ngrid,"zdtsw", 1719 ! . "zdtsw","W.m-2",2, 1720 ! . zdtsw) 1721 1722 1723 !! 1724 !! ***WRF: everything below is kept for reference 1725 !! 1726 ! 1727 !c ========================================================== 1728 !c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing 1729 !c any variable for diagnostic (output with period 1730 !c "ecritphy", set in "run.def") 1731 !c ========================================================== 1732 !c WRITEDIAGFI can ALSO be called from any other subroutines 1733 !c for any variables !! 1734 ! call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, 1735 ! & emis) 1736 ! call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, 1737 ! & tsurf) 1738 ! call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) 1739 ! call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2, 1740 ! & co2ice) 1741 !c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7", 1742 !c & "K",2,zt(1,7)) 1743 ! call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2, 1744 ! & fluxsurf_lw) 1745 ! call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2, 1746 ! & fluxsurf_sw_tot) 1747 ! call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2, 1748 ! & fluxtop_lw) 1749 ! call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2, 1750 ! & fluxtop_sw_tot) 1751 ! call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) 1752 !c call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) 1753 !c call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) 1754 !c call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) 1755 ! call WRITEDIAGFI(ngrid,"rho","density","none",3,rho) 1756 !c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2) 1757 !c call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh) 1758 !c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay) 1759 !c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2, 1760 !c & zstress) 1761 ! 1762 !c ---------------------------------------------------------- 1763 !c Outputs of the CO2 cycle 1764 !c ---------------------------------------------------------- 1765 ! 1766 ! if (tracer.and.(igcm_co2.ne.0)) then 1767 !! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer", 1768 !! & "kg/kg",2,zq(1,1,igcm_co2)) 1769 !! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio", 1770 !! & "kg/kg",3,zq(1,1,igcm_co2)) 1771 ! 1772 ! ! Compute co2 column 1773 ! call zerophys(ngrid,co2col) 1774 ! do l=1,nlayermx 1775 ! do ig=1,ngrid 1776 ! co2col(ig)=co2col(ig)+ 1777 ! & zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g 1778 ! enddo 1779 ! enddo 1780 ! call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2, 1781 ! & co2col) 1782 ! endif ! of if (tracer.and.(igcm_co2.ne.0)) 1783 ! 1784 !c ---------------------------------------------------------- 1785 !c Outputs of the water cycle 1786 !c ---------------------------------------------------------- 1787 ! if (tracer) then 1788 ! if (water) then 1789 ! 1790 ! CALL WRITEDIAGFI(ngridmx,'mtot', 1791 ! & 'total mass of water vapor', 1792 ! & 'kg/m2',2,mtot) 1793 ! CALL WRITEDIAGFI(ngridmx,'icetot', 1794 ! & 'total mass of water ice', 1795 ! & 'kg/m2',2,icetot) 1796 !c If activice is true, tauTES is computed in aeropacity.F; 1797 ! if (.not.activice) then 1798 ! CALL WRITEDIAGFI(ngridmx,'tauTES', 1799 ! & 'tau abs 825 cm-1', 1800 ! & '',2,tauTES) 1801 ! endif 1802 ! 1803 ! call WRITEDIAGFI(ngridmx,'h2o_ice_s', 1804 ! & 'surface h2o_ice', 1805 ! & 'kg.m-2',2,qsurf(1,igcm_h2o_ice)) 1806 ! 1807 ! if (activice) then 1808 !c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate', 1809 !c & 'w.m-2',3,zdtsw) 1810 !c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate', 1811 !c & 'w.m-2',3,zdtlw) 1812 ! endif !(activice) 1813 ! endif !(water) 1814 ! 1815 ! 1816 ! if (water.and..not.photochem) then 1817 ! iq=nq 1818 !c write(str2(1:2),'(i2.2)') iq 1819 !c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud', 1820 !c & 'kg.m-2',2,zdqscloud(1,iq)) 1821 !c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim', 1822 !c & 'kg/kg',3,zdqchim(1,1,iq)) 1823 !c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif', 1824 !c & 'kg/kg',3,zdqdif(1,1,iq)) 1825 !c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj', 1826 !c & 'kg/kg',3,zdqadj(1,1,iq)) 1827 !c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c', 1828 !c & 'kg/kg',3,zdqc(1,1,iq)) 1829 ! endif !(water.and..not.photochem) 1830 ! endif 1831 ! 1832 !c ---------------------------------------------------------- 1833 !c Outputs of the dust cycle 1834 !c ---------------------------------------------------------- 1835 ! 1836 ! call WRITEDIAGFI(ngridmx,'taudustvis', 1837 ! & 'Dust optical depth',' ',2,tau(1,1)) 1838 ! 1839 ! if (tracer.and.(dustbin.ne.0)) then 1840 ! call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1)) 1841 ! if (doubleq) then 1842 ! call WRITEDIAGFI(ngridmx,'qsurf','qsurf', 1843 ! & 'kg.m-2',2,qsurf(1,1)) 1844 ! call WRITEDIAGFI(ngridmx,'Nsurf','N particles', 1845 ! & 'N.m-2',2,qsurf(1,2)) 1846 ! call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift', 1847 ! & 'kg.m-2.s-1',2,zdqsdev(1,1)) 1848 ! call WRITEDIAGFI(ngridmx,'dqssed','sedimentation', 1849 ! & 'kg.m-2.s-1',2,zdqssed(1,1)) 1850 ! do l=1,nlayer 1851 ! do ig=1, ngrid 1852 ! reff(ig,l)= ref_r0 * 1853 ! & (r3n_q*pq(ig,l,1)/max(pq(ig,l,2),0.01))**(1./3.) 1854 ! reff(ig,l)=min(max(reff(ig,l),1.e-10),500.e-6) 1855 ! end do 1856 ! end do 1857 ! call WRITEDIAGFI(ngridmx,'reff','reff','m',3,reff) 1858 ! else 1859 ! do iq=1,dustbin 1860 ! write(str2(1:2),'(i2.2)') iq 1861 ! call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio', 1862 ! & 'kg/kg',3,zq(1,1,iq)) 1863 ! call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf', 1864 ! & 'kg.m-2',2,qsurf(1,iq)) 1865 ! end do 1866 ! endif ! (doubleq) 1867 ! end if ! (tracer.and.(dustbin.ne.0)) 1868 ! 1869 !c ---------------------------------------------------------- 1870 !c Output in netcdf file "diagsoil.nc" for subterranean 1871 !c variables (output every "ecritphy", as for writediagfi) 1872 !c ---------------------------------------------------------- 1873 ! 1874 ! ! Write soil temperature 1875 !! call writediagsoil(ngrid,"soiltemp","Soil temperature","K", 1876 !! & 3,tsoil) 1877 ! ! Write surface temperature 1878 !! call writediagsoil(ngrid,"tsurf","Surface temperature","K", 1879 !! & 2,tsurf) 1880 ! 1881 !c ========================================================== 1882 !c END OF WRITEDIAGFI 1883 !c ========================================================== 1369 #else 1370 1371 c ========================================================== 1372 c WRITEDIAGFI: Outputs in netcdf file "DIAGFI", containing 1373 c any variable for diagnostic (output with period 1374 c "ecritphy", set in "run.def") 1375 c ========================================================== 1376 c WRITEDIAGFI can ALSO be called from any other subroutines 1377 c for any variables !! 1378 c call WRITEDIAGFI(ngrid,"emis","Surface emissivity","w.m-1",2, 1379 c & emis) 1380 ! call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",3,zplay) 1381 ! call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",3,zplev) 1382 call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",2, 1383 & tsurf) 1384 call WRITEDIAGFI(ngrid,"ps","surface pressure","Pa",2,ps) 1385 call WRITEDIAGFI(ngrid,"co2ice","co2 ice thickness","kg.m-2",2, 1386 & co2ice) 1387 c call WRITEDIAGFI(ngrid,"temp7","temperature in layer 7", 1388 c & "K",2,zt(1,7)) 1389 call WRITEDIAGFI(ngrid,"fluxsurf_lw","fluxsurf_lw","W.m-2",2, 1390 & fluxsurf_lw) 1391 call WRITEDIAGFI(ngrid,"fluxsurf_sw","fluxsurf_sw","W.m-2",2, 1392 & fluxsurf_sw_tot) 1393 call WRITEDIAGFI(ngrid,"fluxtop_lw","fluxtop_lw","W.m-2",2, 1394 & fluxtop_lw) 1395 call WRITEDIAGFI(ngrid,"fluxtop_sw","fluxtop_sw","W.m-2",2, 1396 & fluxtop_sw_tot) 1397 call WRITEDIAGFI(ngrid,"temp","temperature","K",3,zt) 1398 call WRITEDIAGFI(ngrid,"u","Zonal wind","m.s-1",3,zu) 1399 call WRITEDIAGFI(ngrid,"v","Meridional wind","m.s-1",3,zv) 1400 call WRITEDIAGFI(ngrid,"w","Vertical wind","m.s-1",3,pw) 1401 call WRITEDIAGFI(ngrid,"rho","density","none",3,rho) 1402 c call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",3,q2) 1403 ! call WRITEDIAGFI(ngrid,'Teta','T potentielle','K',3,zh) 1404 c call WRITEDIAGFI(ngrid,"pressure","Pressure","Pa",3,pplay) 1405 c call WRITEDIAGFI(ngrid,"ssurf","Surface stress","N.m-2",2, 1406 c & zstress) 1407 c call WRITEDIAGFI(ngridmx,'sw_htrt','sw heat. rate', 1408 c & 'w.m-2',3,zdtsw) 1409 c call WRITEDIAGFI(ngridmx,'lw_htrt','lw heat. rate', 1410 c & 'w.m-2',3,zdtlw) 1411 c CALL WRITEDIAGFI(ngridmx,'tauTESap', 1412 c & 'tau abs 825 cm-1', 1413 c & '',2,tauTES) 1414 1415 1416 c ---------------------------------------------------------- 1417 c Outputs of the CO2 cycle 1418 c ---------------------------------------------------------- 1419 1420 if (tracer.and.(igcm_co2.ne.0)) then 1421 ! call WRITEDIAGFI(ngrid,"co2_l1","co2 mix. ratio in 1st layer", 1422 ! & "kg/kg",2,zq(1,1,igcm_co2)) 1423 ! call WRITEDIAGFI(ngrid,"co2","co2 mass mixing ratio", 1424 ! & "kg/kg",3,zq(1,1,igcm_co2)) 1425 1426 ! Compute co2 column 1427 call zerophys(ngrid,co2col) 1428 do l=1,nlayermx 1429 do ig=1,ngrid 1430 co2col(ig)=co2col(ig)+ 1431 & zq(ig,l,igcm_co2)*(pplev(ig,l)-pplev(ig,l+1))/g 1432 enddo 1433 enddo 1434 call WRITEDIAGFI(ngrid,"co2col","CO2 column","kg.m-2",2, 1435 & co2col) 1436 endif ! of if (tracer.and.(igcm_co2.ne.0)) 1437 1438 c ---------------------------------------------------------- 1439 c Outputs of the water cycle 1440 c ---------------------------------------------------------- 1441 if (tracer) then 1442 if (water) then 1443 1444 CALL WRITEDIAGFI(ngridmx,'mtot', 1445 & 'total mass of water vapor', 1446 & 'kg/m2',2,mtot) 1447 CALL WRITEDIAGFI(ngridmx,'icetot', 1448 & 'total mass of water ice', 1449 & 'kg/m2',2,icetot) 1450 c vmr=zq(1:ngridmx,1:nlayermx,igcm_h2o_ice) 1451 c & *mugaz/mmol(igcm_h2o_ice) 1452 c call WRITEDIAGFI(ngridmx,'vmr_h2oice','h2o ice vmr', 1453 c & 'mol/mol',3,vmr) 1454 CALL WRITEDIAGFI(ngridmx,'reffice', 1455 & 'Mean reff', 1456 & 'm',2,rave) 1457 c call WRITEDIAGFI(ngridmx,'rice','Ice particle size', 1458 c & 'm',3,rice) 1459 c If activice is true, tauTES is computed in aeropacity.F; 1460 if (.not.activice) then 1461 CALL WRITEDIAGFI(ngridmx,'tauTESap', 1462 & 'tau abs 825 cm-1', 1463 & '',2,tauTES) 1464 endif 1465 call WRITEDIAGFI(ngridmx,'h2o_ice_s', 1466 & 'surface h2o_ice', 1467 & 'kg.m-2',2,qsurf(1,igcm_h2o_ice)) 1468 endif !(water) 1469 1470 1471 if (water.and..not.photochem) then 1472 iq=nq 1473 c write(str2(1:2),'(i2.2)') iq 1474 c call WRITEDIAGFI(ngridmx,'dqs'//str2,'dqscloud', 1475 c & 'kg.m-2',2,zdqscloud(1,iq)) 1476 c call WRITEDIAGFI(ngridmx,'dqch'//str2,'var chim', 1477 c & 'kg/kg',3,zdqchim(1,1,iq)) 1478 c call WRITEDIAGFI(ngridmx,'dqd'//str2,'var dif', 1479 c & 'kg/kg',3,zdqdif(1,1,iq)) 1480 c call WRITEDIAGFI(ngridmx,'dqa'//str2,'var adj', 1481 c & 'kg/kg',3,zdqadj(1,1,iq)) 1482 c call WRITEDIAGFI(ngridmx,'dqc'//str2,'var c', 1483 c & 'kg/kg',3,zdqc(1,1,iq)) 1484 endif !(water.and..not.photochem) 1485 endif 1486 1487 c ---------------------------------------------------------- 1488 c Outputs of the dust cycle 1489 c ---------------------------------------------------------- 1490 1491 c call WRITEDIAGFI(ngridmx,'tauref', 1492 c & 'Dust ref opt depth','NU',2,tauref) 1493 1494 if (tracer.and.(dustbin.ne.0)) then 1495 c call WRITEDIAGFI(ngridmx,'tau','taudust','SI',2,tau(1,1)) 1496 if (doubleq) then 1497 c call WRITEDIAGFI(ngridmx,'qsurf','qsurf', 1498 c & 'kg.m-2',2,qsurf(1,1)) 1499 c call WRITEDIAGFI(ngridmx,'Nsurf','N particles', 1500 c & 'N.m-2',2,qsurf(1,2)) 1501 c call WRITEDIAGFI(ngridmx,'dqsdev','ddevil lift', 1502 c & 'kg.m-2.s-1',2,zdqsdev(1,1)) 1503 c call WRITEDIAGFI(ngridmx,'dqssed','sedimentation', 1504 c & 'kg.m-2.s-1',2,zdqssed(1,1)) 1505 call WRITEDIAGFI(ngridmx,'reffdust','reffdust', 1506 & 'm',3,rdust*ref_r0) 1507 call WRITEDIAGFI(ngridmx,'dustq','Dust mass mr', 1508 & 'kg/kg',3,pq(1,1,igcm_dust_mass)) 1509 c call WRITEDIAGFI(ngridmx,'dustN','Dust number', 1510 c & 'part/kg',3,pq(1,1,igcm_dust_number)) 1511 else 1512 do iq=1,dustbin 1513 write(str2(1:2),'(i2.2)') iq 1514 call WRITEDIAGFI(ngridmx,'q'//str2,'mix. ratio', 1515 & 'kg/kg',3,zq(1,1,iq)) 1516 call WRITEDIAGFI(ngridmx,'qsurf'//str2,'qsurf', 1517 & 'kg.m-2',2,qsurf(1,iq)) 1518 end do 1519 endif ! (doubleq) 1520 c if (submicron) then 1521 c call WRITEDIAGFI(ngridmx,'dustsubm','subm mass mr', 1522 c & 'kg/kg',3,pq(1,1,igcm_dust_submicron)) 1523 c endif ! (submicron) 1524 end if ! (tracer.and.(dustbin.ne.0)) 1525 1526 c ---------------------------------------------------------- 1527 c Outputs of thermals 1528 c ---------------------------------------------------------- 1529 if (calltherm) then 1530 1531 ! call WRITEDIAGFI(ngrid,'dtke', 1532 ! & 'tendance tke thermiques','m**2/s**2', 1533 ! & 3,dtke_th) 1534 ! call WRITEDIAGFI(ngrid,'d_u_ajs', 1535 ! & 'tendance u thermiques','m/s', 1536 ! & 3,pdu_th*ptimestep) 1537 ! call WRITEDIAGFI(ngrid,'d_v_ajs', 1538 ! & 'tendance v thermiques','m/s', 1539 ! & 3,pdv_th*ptimestep) 1540 ! if (tracer) then 1541 ! if (nq .eq. 2) then 1542 ! call WRITEDIAGFI(ngrid,'deltaq_th', 1543 ! & 'delta q thermiques','kg/kg', 1544 ! & 3,ptimestep*pdq_th(:,:,2)) 1545 ! endif 1546 ! endif 1547 1548 lmax_th_out(:)=real(lmax_th(:)) 1549 1550 call WRITEDIAGFI(ngridmx,'lmax_th', 1551 & 'hauteur du thermique','K', 1552 & 2,lmax_th_out) 1553 call WRITEDIAGFI(ngridmx,'hfmax_th', 1554 & 'maximum TH heat flux','K.m/s', 1555 & 2,hfmax_th) 1556 call WRITEDIAGFI(ngridmx,'wmax_th', 1557 & 'maximum TH vertical velocity','m/s', 1558 & 2,wmax_th) 1559 1560 endif 1561 1562 c ---------------------------------------------------------- 1563 c Output in netcdf file "diagsoil.nc" for subterranean 1564 c variables (output every "ecritphy", as for writediagfi) 1565 c ---------------------------------------------------------- 1566 1567 ! Write soil temperature 1568 ! call writediagsoil(ngrid,"soiltemp","Soil temperature","K", 1569 ! & 3,tsoil) 1570 ! Write surface temperature 1571 ! call writediagsoil(ngrid,"tsurf","Surface temperature","K", 1572 ! & 2,tsurf) 1573 1574 c ========================================================== 1575 c END OF WRITEDIAGFI 1576 c ========================================================== 1577 #endif 1884 1578 1885 1579 ELSE ! if(ngrid.eq.1) 1886 1580 1581 #ifndef MESOSCALE 1887 1582 print*,'Ls =',zls*180./pi, 1888 1583 & ' tauref(700 Pa) =',tauref … … 1901 1596 c CALL writeg1d(ngrid,nlayer,pw,'w','m.s-1') 1902 1597 1903 !! or output in diagfi.nc (for testphys1d) 1904 ! call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps) 1905 ! call WRITEDIAGFI(ngridmx,'temp','Temperature', 1906 ! & 'K',1,zt) 1907 ! 1908 ! if(tracer) then 1909 !c CALL writeg1d(ngrid,1,tau,'tau','SI') 1910 ! do iq=1,nq 1911 !c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 1912 ! call WRITEDIAGFI(ngridmx,trim(noms(iq)), 1913 ! & trim(noms(iq)),'kg/kg',1,zq(1,1,iq)) 1914 ! end do 1915 ! end if 1916 ! 1917 ! zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g 1918 ! 1919 ! do l=2,nlayer-1 1920 ! tmean=zt(1,l) 1921 ! if(zt(1,l).ne.zt(1,l-1)) 1922 ! & tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1)) 1923 ! zlocal(l)= zlocal(l-1) 1924 ! & -log(pplay(1,l)/pplay(1,l-1))*rnew(1,l)*tmean/g 1925 ! enddo 1926 ! zlocal(nlayer)= zlocal(nlayer-1)- 1927 ! & log(pplay(1,nlayer)/pplay(1,nlayer-1))* 1928 ! & rnew(1,nlayer)*tmean/g 1598 ! THERMALS STUFF 1D 1599 1600 if(calltherm) then 1601 1602 lmax_th_out(:)=real(lmax_th(:)) 1603 1604 call WRITEDIAGFI(ngridmx,'lmax_th', 1605 & 'hauteur du thermique','point', 1606 & 0,lmax_th_out) 1607 call WRITEDIAGFI(ngridmx,'hfmax_th', 1608 & 'maximum TH heat flux','K.m/s', 1609 & 0,hfmax_th) 1610 call WRITEDIAGFI(ngridmx,'wmax_th', 1611 & 'maximum TH vertical velocity','m/s', 1612 & 0,wmax_th) 1613 1614 1615 co2col(:)=0. 1616 dummycol(:)=0. 1617 if (tracer) then 1618 do l=1,nlayermx 1619 do ig=1,ngrid 1620 co2col(ig)=co2col(ig)+ 1621 & zq(ig,l,1)*(pplev(ig,l)-pplev(ig,l+1))/g 1622 if (nqmx .gt. 1) then 1623 iq=2 ! to avoid out-of-bounds spotting by picky compilers 1624 ! when gcm is compiled with only one tracer 1625 dummycol(ig)=dummycol(ig)+ 1626 & zq(ig,l,iq)*(pplev(ig,l)-pplev(ig,l+1))/g 1627 endif 1628 enddo 1629 enddo 1630 1631 end if 1632 call WRITEDIAGFI(ngrid,'co2col','integrated co2 mass' & 1633 & ,'kg/m-2',0,co2col) 1634 call WRITEDIAGFI(ngrid,'dummycol','integrated dummy mass' & 1635 & ,'kg/m-2',0,dummycol) 1636 endif 1637 call WRITEDIAGFI(ngrid,'w','vertical velocity' & 1638 & ,'m/s',1,pw) 1639 call WRITEDIAGFI(ngrid,"q2","q2","kg.m-3",1,q2) 1640 call WRITEDIAGFI(ngrid,"tsurf","Surface temperature","K",0, 1641 & tsurf) 1642 1643 call WRITEDIAGFI(ngrid,"pplay","Pressure","Pa",1,zplay) 1644 call WRITEDIAGFI(ngrid,"pplev","Pressure","Pa",1,zplev) 1645 ! or output in diagfi.nc (for testphys1d) 1646 call WRITEDIAGFI(ngridmx,'ps','Surface pressure','Pa',0,ps) 1647 call WRITEDIAGFI(ngridmx,'temp','Temperature', 1648 & 'K',1,zt) 1649 1650 if(tracer) then 1651 c CALL writeg1d(ngrid,1,tau,'tau','SI') 1652 do iq=1,nq 1653 c CALL writeg1d(ngrid,nlayer,zq(1,1,iq),noms(iq),'kg/kg') 1654 call WRITEDIAGFI(ngridmx,trim(noms(iq)), 1655 & trim(noms(iq)),'kg/kg',1,zq(1,1,iq)) 1656 end do 1657 end if 1658 1659 zlocal(1)=-log(pplay(1,1)/pplev(1,1))* Rnew(1,1)*zt(1,1)/g 1660 1661 do l=2,nlayer-1 1662 tmean=zt(1,l) 1663 if(zt(1,l).ne.zt(1,l-1)) 1664 & tmean=(zt(1,l)-zt(1,l-1))/log(zt(1,l)/zt(1,l-1)) 1665 zlocal(l)= zlocal(l-1) 1666 & -log(pplay(1,l)/pplay(1,l-1))*rnew(1,l)*tmean/g 1667 enddo 1668 zlocal(nlayer)= zlocal(nlayer-1)- 1669 & log(pplay(1,nlayer)/pplay(1,nlayer-1))* 1670 & rnew(1,nlayer)*tmean/g 1671 #endif 1929 1672 1930 1673 END IF ! if(ngrid.ne.1) 1931 1674 1932 1675 icount=icount+1 1676 #ifdef MESOSCALE 1933 1677 write(*,*) 'now, back to the dynamical core...' 1934 1678 #endif -
trunk/LMDZ.MARS/libf/phymars/physiq.F
r224 r226 1 SUBROUTINE physiq(ngrid,nlayer,nq, 2 $ firstcall,lastcall, 3 $ pday,ptime,ptimestep, 4 $ pplev,pplay,pphi, 5 $ pu,pv,pt,pq, 6 $ pw, 7 $ pdu,pdv,pdt,pdq,pdpsrf,tracerdyn) 1 SUBROUTINE physiq( 2 $ ngrid,nlayer,nq 3 $ ,firstcall,lastcall 4 $ ,pday,ptime,ptimestep 5 $ ,pplev,pplay,pphi 6 $ ,pu,pv,pt,pq 7 $ ,pw 8 $ ,pdu,pdv,pdt,pdq,pdpsrf,tracerdyn 9 $ ) 8 10 9 11 IMPLICIT NONE … … 206 208 207 209 CHARACTER*80 fichier 208 INTEGER l,ig,ierr,igout,iq,i, tapphys,k 209 LOGICAL end_of_file 210 INTEGER l,ig,ierr,igout,iq,i, tapphys 210 211 211 212 REAL fluxsurf_lw(ngridmx) !incident LW (IR) surface flux (W.m-2) … … 1080 1081 c Sum of fluxes in solar spectral bands (for output only) 1081 1082 DO ig=1,ngrid 1082 1083 1083 fluxtop_sw_tot(ig)=fluxtop_sw(ig,1) + fluxtop_sw(ig,2) 1084 fluxsurf_sw_tot(ig)=fluxsurf_sw(ig,1) + fluxsurf_sw(ig,2) 1084 1085 ENDDO 1085 1086 c ******* TEST ****************************************************** … … 1121 1122 1122 1123 IF (ngrid.NE.1) THEN 1123 print*,'Ls =',zls*180./pi ,1124 & ' tauref(700 Pa,lat=0) =',tauref(ngrid/2),1125 & ' tau(Viking1) =',tau(ig_vl1,1)1124 print*,'Ls =',zls*180./pi 1125 & ,' tauref(700 Pa,lat=0) =',tauref(ngrid/2) 1126 & ,' tau(Viking1) =',tau(ig_vl1,1) 1126 1127 1127 1128 -
trunk/LMDZ.MARS/libf/phymars/surfdat.h
r224 r226 19 19 real iceradius , dtemisice 20 20 real zmea,zstd,zsig,zgam,zthe 21 real z0 ! surface roughness leng ht(m)21 real z0 ! surface roughness length (m) 22 22 real z0_default ! default (constant over planet) surface roughness (m)
Note: See TracChangeset
for help on using the changeset viewer.
