source: lmdz_wrf/trunk/WRFV3/lmdz/phyaqua_mod.F90

Last change on this file was 1, checked in by lfita, 10 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 28.2 KB
Line 
1! Routines complementaires pour la physique planetaire.
2
3MODULE phyaqua_mod
4
5  CONTAINS
6
7      subroutine iniaqua(nlon,latfi,lonfi,iflag_phys)
8
9!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
10!  Creation d'un etat initial et de conditions aux limites
11!  (resp startphy.nc et limit.nc) pour des configurations idealisees
12! du modele LMDZ dans sa version terrestre.
13!  iflag_phys est un parametre qui controle
14!  iflag_phys = N 
15!    de 100 a 199 : aqua planetes avec SST forcees
16!                 N-100 determine le type de SSTs
17!    de 200 a 299 : terra planetes avec Ts calcule
18!       
19!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
20
21      USE comgeomphy, only : rlatd,rlond
22      USE dimphy, only : klon
23      USE surface_data, only : type_ocean,ok_veget
24      USE pbl_surface_mod, only : pbl_surface_init
25      USE fonte_neige_mod, only : fonte_neige_init
26      USE phys_state_var_mod
27      USE control_mod, only : dayref,nday,iphysiq
28      USE indice_sol_mod
29
30      USE IOIPSL
31      IMPLICIT NONE
32
33#include "dimensions.h"
34!   #include "dimphy.h"
35!   #include "YOMCST.h"
36#include "comconst.h"
37#include "clesphys.h"
38#include "dimsoil.h"
39#include "temps.h"
40
41      integer,intent(in) :: nlon,iflag_phys
42!IM ajout latfi, lonfi
43      real,intent(in) :: lonfi(nlon),latfi(nlon)
44
45      INTEGER type_profil,type_aqua
46
47!c  Ajouts initialisation des surfaces
48      REAL :: run_off_lic_0(nlon)
49      REAL :: qsolsrf(nlon,nbsrf),snsrf(nlon,nbsrf)
50      REAL :: frugs(nlon,nbsrf)
51      REAL :: agesno(nlon,nbsrf)
52      REAL :: tsoil(nlon,nsoilmx,nbsrf)
53      REAL :: tslab(nlon), seaice(nlon)
54      REAL evap(nlon,nbsrf),fder(nlon)
55
56
57
58!c    Arguments :
59!c    -----------
60
61!      integer radpas 
62      integer it,unit,i,k,itap
63
64      real airefi,zcufi,zcvfi
65
66      real rugos,albedo
67      REAL tsurf
68      REAL time,timestep,day,day0
69      real qsol_f,qsol(nlon)
70      real rugsrel(nlon)
71!      real zmea(nlon),zstd(nlon),zsig(nlon)
72!      real zgam(nlon),zthe(nlon),zpic(nlon),zval(nlon)
73!      real rlon(nlon),rlat(nlon)
74      logical alb_ocean
75!      integer demih_pas
76
77      CHARACTER*80 ans,file_forctl, file_fordat, file_start
78      character*100 file,var
79      character*2 cnbl
80
81      REAL phy_nat(nlon,360)
82      REAL phy_alb(nlon,360)
83      REAL phy_sst(nlon,360)
84      REAL phy_bil(nlon,360)
85      REAL phy_rug(nlon,360)
86      REAL phy_ice(nlon,360)
87      REAL phy_fter(nlon,360)
88      REAL phy_foce(nlon,360)
89      REAL phy_fsic(nlon,360)
90      REAL phy_flic(nlon,360)
91
92      integer, save::  read_climoz=0 ! read ozone climatology
93
94! intermediate variables to use getin
95      integer :: nbapp_rad_omp
96      real :: co2_ppm_omp,solaire_omp
97      logical :: alb_ocean_omp
98      real :: rugos_omp
99!-------------------------------------------------------------------------
100!  declaration pour l'appel a phyredem
101!-------------------------------------------------------------------------
102
103!      real pctsrf(nlon,nbsrf),ftsol(nlon,nbsrf)
104      real falbe(nlon,nbsrf),falblw(nlon,nbsrf)
105!      real pbl_tke(nlon,llm,nbsrf)
106!      real rain_fall(nlon),snow_fall(nlon)
107!      real solsw(nlon), sollw(nlon),radsol(nlon)
108!      real t_ancien(nlon,llm),q_ancien(nlon,llm),rnebcon(nlon,llm)
109!      real ratqs(nlon,llm)
110!      real clwcon(nlon,llm)
111
112      INTEGER        longcles
113      PARAMETER    ( longcles = 20 )
114      REAL clesphy0( longcles )
115
116
117!c-----------------------------------------------------------------------
118!c   dynamial tendencies :
119!c   ---------------------
120
121      INTEGER l,ierr,aslun
122
123      REAL longitude,latitude
124      REAL paire
125
126      DATA latitude,longitude/48.,0./
127
128!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
129! INITIALISATIONS
130!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
131
132!-----------------------------------------------------------------------
133!    Initialisations  des constantes
134!    -------------------------------
135
136
137      type_aqua=iflag_phys/100
138      type_profil=iflag_phys-type_aqua*100
139      print*,'iniaqua:type_aqua, type_profil',type_aqua, type_profil
140
141      if (klon.ne.nlon) then
142        write(*,*)"iniaqua: klon=",klon," nlon=",nlon
143        stop'probleme de dimensions dans iniaqua'
144      endif
145      call phys_state_var_init(read_climoz)
146
147
148      read_climoz=0
149      day0=217.
150      day=day0
151      it=0
152      time=0.
153
154!IM ajout latfi, lonfi
155      rlatd=latfi
156      rlond=lonfi
157      rlat=rlatd*180./pi
158      rlon=rlond*180./pi
159
160!-----------------------------------------------------------------------
161!  initialisations de la physique
162!-----------------------------------------------------------------------
163
164         day_ini=dayref
165         day_end=day_ini+nday
166         airefi=1.
167         zcufi=1.
168         zcvfi=1.
169!$OMP MASTER
170      nbapp_rad_omp=24
171      CALL getin('nbapp_rad',nbapp_rad_omp)
172!$OMP END MASTER
173!$OMP BARRIER
174      nbapp_rad=nbapp_rad_omp
175
176!---------------------------------------------------------------------
177!c Creation des conditions aux limites:
178!c ------------------------------------
179! Initialisations des constantes
180! Ajouter les manquants dans planete.def... (albedo etc)
181!$OMP MASTER
182      co2_ppm_omp=348.
183      CALL getin('co2_ppm',co2_ppm_omp)
184      solaire_omp=1365.
185      CALL getin('solaire',solaire_omp)
186!      CALL getin('albedo',albedo) ! albedo is set below, depending on type_aqua
187      alb_ocean_omp=.true.
188      CALL getin('alb_ocean',alb_ocean_omp)
189!$OMP END MASTER
190!$OMP BARRIER
191      co2_ppm=co2_ppm_omp
192      solaire=solaire_omp
193      alb_ocean=alb_ocean_omp
194
195      radsol=0.
196      qsol_f=10.
197
198!c  Conditions aux limites:
199!c  -----------------------
200
201      qsol(:)    = qsol_f
202      rugsrel = 0.0    ! (rugsrel = rugoro)
203      rugoro = 0.0
204      u_ancien = 0.0
205      v_ancien = 0.0
206      agesno  = 50.0
207! Relief plat
208      zmea = 0.
209      zstd = 0.
210      zsig = 0.
211      zgam = 0.
212      zthe = 0.
213      zpic = 0.
214      zval = 0.
215
216! Une seule surface
217      pctsrf=0.
218      if (type_aqua==1) then
219         rugos=1.e-4
220         albedo=0.19
221         pctsrf(:,is_oce)=1.
222      else if (type_aqua==2) then
223         rugos=0.03
224         albedo=0.1
225         pctsrf(:,is_ter)=1.
226      endif
227
228!$OMP MASTER
229      rugos_omp=rugos
230      CALL getin('rugos',rugos_omp)
231!$OMP END MASTER
232!$OMP BARRIER
233      rugos=rugos_omp
234      zmasq(:)=pctsrf(:,is_oce)
235
236! pctsrf_pot(:,is_oce) = 1. - zmasq(:)
237! pctsrf_pot(:,is_sic) = 1. - zmasq(:)
238
239! Si alb_ocean on calcule un albedo oceanique moyen
240!  if (alb_ocean) then
241! Voir pourquoi on avait ca.
242!          CALL ini_alb_oce(phy_alb)
243!      else
244      phy_alb(:,:) = albedo ! albedo land only (old value condsurf_jyg=0.3)
245!      endif !alb_ocean
246     
247      do i=1,360
248!IM Terraplanete   phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
249!IM ajout calcul profil sst selon le cas considere (cf. FBr)
250
251      phy_nat(:,i) = 1.0    ! 0=ocean libre, 1=land, 2=glacier, 3=banquise
252      phy_bil(:,i) = 1.0    ! ne sert que pour les slab_ocean
253      phy_rug(:,i) = rugos  ! longueur rugosite utilisee sur land only
254      phy_ice(:,i) = 0.0    ! fraction de glace (?)
255      phy_fter(:,i) = pctsrf(:,is_ter)  ! fraction de glace (?)
256      phy_foce(:,i) = pctsrf(:,is_oce)  ! fraction de glace (?)
257      phy_fsic(:,i) = pctsrf(:,is_sic)  ! fraction de glace (?)
258      phy_flic(:,i) = pctsrf(:,is_lic)  ! fraction de glace (?)
259      enddo
260!IM calcul profil sst
261      call profil_sst(nlon, rlatd, type_profil, phy_sst)
262
263      call writelim                                                                  &
264       &   (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,                    &
265       &    phy_fter,phy_foce,phy_flic,phy_fsic)
266
267
268!---------------------------------------------------------------------
269!c Ecriture de l'etat initial:
270!c ---------------------------
271
272!C
273!C Ecriture etat initial physique
274!C
275      timestep   = dtvr * FLOAT(iphysiq)
276      radpas    = NINT (daysec/timestep/ FLOAT(nbapp_rad) )
277
278      DO i = 1, longcles
279       clesphy0(i) = 0.
280      ENDDO
281      clesphy0(1) = FLOAT( iflag_con )
282      clesphy0(2) = FLOAT( nbapp_rad )
283!c     IF( cycle_diurne  ) clesphy0(3) =  1.
284      clesphy0(3)=1. ! cycle_diurne
285      clesphy0(4)=1. ! soil_model
286      clesphy0(5)=1. ! new_oliq
287      clesphy0(6)=0. ! ok_orodr
288      clesphy0(7)=0. ! ok_orolf
289      clesphy0(8)=0. ! ok_limitvrai
290
291
292!c=======================================================================
293!c  Profils initiaux
294!c=======================================================================
295
296! On initialise les temperatures de surfaces comme les sst
297      do i=1,nlon
298         ftsol(i,:)=phy_sst(i,1)
299         tsoil(i,:,:)=phy_sst(i,1)
300         tslab(i)=phy_sst(i,1)
301      enddo
302
303      falbe(:,:)=albedo
304      falblw(:,:)=albedo
305      rain_fall(:)=0.
306      snow_fall(:)=0.
307      solsw(:)=0.
308      sollw(:)=0.
309      radsol(:)=0.
310
311!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
312!  intialisation bidon mais pas grave
313      t_ancien(:,:)=0.
314      q_ancien(:,:)=0.
315!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
316      rnebcon=0.
317      ratqs=0.
318      clwcon=0.
319      pbl_tke=1.e-8
320
321! variables supplementaires pour appel a plb_surface_init
322      fder(:)=0.
323      seaice(:)=0.
324      run_off_lic_0=0.
325      evap=0.
326
327
328! Initialisations necessaires avant phyredem
329      type_ocean = "force"
330      call fonte_neige_init(run_off_lic_0)
331      qsolsrf(:,:)=qsol(1) ! humidite du sol des sous surface
332      snsrf(:,:)=0.        ! couverture de neige des sous surface
333      frugs(:,:)=rugos        ! couverture de neige des sous surface
334
335
336      call pbl_surface_init(qsol, fder, snsrf, qsolsrf,                              &
337       &     evap, frugs, agesno, tsoil)
338
339        print*,'iniaqua: before phyredem'
340
341      falb1=albedo
342      falb2=albedo
343      zmax0=0.
344      f0=0.
345      ema_work1=0.
346      ema_work2=0.
347      wake_deltat=0.
348      wake_deltaq=0.
349      wake_s=0.
350      wake_cstar=0.
351      wake_pe=0.
352      wake_fip=0.
353      fm_therm=0.
354      entr_therm=0.
355      detr_therm=0.
356
357
358      CALL phyredem ("startphy.nc")
359
360        print*,'iniaqua: after phyredem'
361      call phys_state_var_end
362
363      return
364
365      end SUBROUTINE iniaqua
366
367!c====================================================================
368!c====================================================================
369      SUBROUTINE zenang_an(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
370      USE dimphy
371      IMPLICIT none
372!c====================================================================
373!c=============================================================
374!c         CALL zenang(cycle_diurne,gmtime,rlat,rlon,rmu0,fract)
375!c Auteur : A. Campoy et F. Hourdin
376!c Objet  : calculer les valeurs moyennes du cos de l'angle zenithal
377!c          et l'ensoleillement moyen entre gmtime1 et gmtime2
378!c          connaissant la declinaison, la latitude et la longitude.
379!c
380!c   Dans cette version particuliere, on calcule le rayonnement
381!c  moyen sur l'année à chaque latitude.
382!c   angle zenithal calculé pour obtenir un
383!c   Fit polynomial de  l'ensoleillement moyen au sommet de l'atmosphere
384!c   en moyenne annuelle.
385!c   Spécifique de la terre. Utilisé pour les aqua planetes.
386!c
387!c Rque   : Different de la routine angle en ce sens que zenang
388!c          fournit des moyennes de pmu0 et non des valeurs
389!c          instantanees, du coup frac prend toutes les valeurs
390!c          entre 0 et 1.
391!c Date   : premiere version le 13 decembre 1994
392!c          revu pour  GCM  le 30 septembre 1996
393!c===============================================================
394!c longi----INPUT : la longitude vraie de la terre dans son plan
395!c                  solaire a partir de l'equinoxe de printemps (degre)
396!c gmtime---INPUT : temps universel en fraction de jour
397!c pdtrad---INPUT : pas de temps du rayonnement (secondes)
398!c lat------INPUT : latitude en degres
399!c long-----INPUT : longitude en degres
400!c pmu0-----OUTPUT: angle zenithal moyen entre gmtime et gmtime+pdtrad
401!c frac-----OUTPUT: ensoleillement moyen entre gmtime et gmtime+pdtrad
402!c================================================================
403#include "YOMCST.h"
404!c================================================================
405      logical cycle_diurne
406      real  gmtime
407      real rlat(klon), rlon(klon), rmu0(klon), fract(klon)
408!c================================================================
409      integer i
410      real gmtime1, gmtime2
411      real pi_local
412
413
414      real rmu0m(klon),rmu0a(klon)
415!c
416
417      pi_local = 4.0 * ATAN(1.0)
418
419!c================================================================
420!c  Calcul de l'angle zenithal moyen sur la journee
421!c================================================================
422
423      DO i=1,klon
424        fract(i)=1.
425!  Calcule du flux moyen
426        IF (abs(rlat(i)).LE.28.75) THEN
427           rmu0m(i)=(210.1924+206.6059*cos(0.0174533*rlat(i))**2)/1365.
428        ELSEIF (abs(rlat(i)).LE.43.75) THEN
429          rmu0m(i)=(187.4562+236.1853*cos(0.0174533*rlat(i))**2)/1365.
430        ELSEIF (abs(rlat(i)).LE.71.25) THEN
431          rmu0m(i)=(162.4439+284.1192*cos(0.0174533*rlat(i))**2)/1365.
432        ELSE
433          rmu0m(i)=(172.8125+183.7673*cos(0.0174533*rlat(i))**2)/1365.
434        ENDIF
435      ENDDO
436
437!c================================================================
438!  Avec ou sans cycle diurne
439!c================================================================
440
441      IF (cycle_diurne) THEN
442
443!  On redecompose flux  au sommet suivant un cycle diurne idealise
444!  identique a toutes les latitudes.
445
446         DO i=1,klon
447           rmu0a(i)=2.*rmu0m(i)*sqrt(2.)*pi_local/(4.-pi_local)
448           rmu0(i)=rmu0a(i)*abs(sin(pi_local*gmtime+pi_local*                        &
449       &      rlon(i)/360.))-rmu0a(i)/sqrt(2.)
450         ENDDO
451
452         DO i=1,klon
453           IF (rmu0(i).LE.0.) THEN
454              rmu0(i)=0.
455              fract(i)=0.
456           ELSE
457              fract(i)=1.
458           ENDIF
459         ENDDO
460
461!  Affichage de l'angel zenitale
462!     print*,'************************************'
463!     print*,'************************************'
464!     print*,'************************************'
465!     print*,'latitude=',rlat(i),'longitude=',rlon(i)
466!     print*,'rmu0m=',rmu0m(i)
467!     print*,'rmu0a=',rmu0a(i)
468!     print*,'rmu0=',rmu0(i)
469                                                               
470      ELSE
471
472        DO i=1,klon
473           fract(i)=0.5
474           rmu0(i)=rmu0m(i)*2.
475        ENDDO
476
477      ENDIF
478
479      RETURN
480      END SUBROUTINE zenang_an
481
482!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
483
484      subroutine writelim                                                            &
485       &   (klon,phy_nat,phy_alb,phy_sst,phy_bil,phy_rug,phy_ice,                    &
486       &    phy_fter,phy_foce,phy_flic,phy_fsic)
487!c
488      use mod_phys_lmdz_para, only: is_mpi_root,is_omp_root
489      use mod_grid_phy_lmdz, only : klon_glo
490      use mod_phys_lmdz_transfert_para, only : gather
491
492!L. Fita July 2013
493      use netcdf
494!#include "dimensions.h"
495!#include "dimphy.h"
496 
497      integer,intent(in) :: klon
498      real,intent(in) :: phy_nat(klon,360)
499      real,intent(in) :: phy_alb(klon,360)
500      real,intent(in) :: phy_sst(klon,360)
501      real,intent(in) :: phy_bil(klon,360)
502      real,intent(in) :: phy_rug(klon,360)
503      real,intent(in) :: phy_ice(klon,360)
504      real,intent(in) :: phy_fter(klon,360)
505      real,intent(in) :: phy_foce(klon,360)
506      real,intent(in) :: phy_flic(klon,360)
507      real,intent(in) :: phy_fsic(klon,360)
508
509      real :: phy_glo(klon_glo,360) ! temporary variable, to store phy_***(:)
510                                    ! on the whole physics grid
511      INTEGER ierr
512      INTEGER dimfirst(3)
513      INTEGER dimlast(3)
514!c
515      INTEGER nid, ndim, ntim
516      INTEGER dims(2), debut(2), epais(2)
517      INTEGER id_tim
518      INTEGER id_NAT, id_SST, id_BILS, id_RUG, id_ALB
519      INTEGER id_FTER,id_FOCE,id_FSIC,id_FLIC
520 
521      if (is_mpi_root.and.is_omp_root) then
522     
523        PRINT*, 'writelim: Ecriture du fichier limit'
524!c
525        ierr = NF_CREATE ("limit.nc", NF_CLOBBER, nid)
526!c
527        ierr = NF_PUT_ATT_TEXT (nid, NF_GLOBAL, "title", 30,                         &
528       &                       "Fichier conditions aux limites")
529!!        ierr = NF_DEF_DIM (nid, "points_physiques", klon, ndim)
530        ierr = NF_DEF_DIM (nid, "points_physiques", klon_glo, ndim)
531        ierr = NF_DEF_DIM (nid, "time", NF_UNLIMITED, ntim)
532!c
533        dims(1) = ndim
534        dims(2) = ntim
535!c
536!ccc      ierr = NF_DEF_VAR (nid, "TEMPS", NF_DOUBLE, 1,ntim, id_tim)
537        ierr = NF_DEF_VAR (nid, "TEMPS", NF_FLOAT, 1,ntim, id_tim)
538        ierr = NF_PUT_ATT_TEXT (nid, id_tim, "title", 17,                            &
539       &                        "Jour dans l annee")
540!ccc      ierr = NF_DEF_VAR (nid, "NAT", NF_DOUBLE, 2,dims, id_NAT)
541        ierr = NF_DEF_VAR (nid, "NAT", NF_FLOAT, 2,dims, id_NAT)
542        ierr = NF_PUT_ATT_TEXT (nid, id_NAT, "title", 23,                            &
543       &                        "Nature du sol (0,1,2,3)")
544!ccc      ierr = NF_DEF_VAR (nid, "SST", NF_DOUBLE, 2,dims, id_SST)
545        ierr = NF_DEF_VAR (nid, "SST", NF_FLOAT, 2,dims, id_SST)
546        ierr = NF_PUT_ATT_TEXT (nid, id_SST, "title", 35,                            &
547       &                        "Temperature superficielle de la mer")
548!ccc      ierr = NF_DEF_VAR (nid, "BILS", NF_DOUBLE, 2,dims, id_BILS)
549        ierr = NF_DEF_VAR (nid, "BILS", NF_FLOAT, 2,dims, id_BILS)
550        ierr = NF_PUT_ATT_TEXT (nid, id_BILS, "title", 32,                           &
551       &                        "Reference flux de chaleur au sol")
552!ccc      ierr = NF_DEF_VAR (nid, "ALB", NF_DOUBLE, 2,dims, id_ALB)
553        ierr = NF_DEF_VAR (nid, "ALB", NF_FLOAT, 2,dims, id_ALB)
554        ierr = NF_PUT_ATT_TEXT (nid, id_ALB, "title", 19,                            &
555       &                        "Albedo a la surface")
556!ccc      ierr = NF_DEF_VAR (nid, "RUG", NF_DOUBLE, 2,dims, id_RUG)
557        ierr = NF_DEF_VAR (nid, "RUG", NF_FLOAT, 2,dims, id_RUG)
558        ierr = NF_PUT_ATT_TEXT (nid, id_RUG, "title", 8,                             &
559       &                        "Rugosite")
560
561        ierr = NF_DEF_VAR (nid, "FTER", NF_FLOAT, 2,dims, id_FTER)
562        ierr = NF_PUT_ATT_TEXT (nid, id_FTER, "title", 8,"Frac. Terre")
563        ierr = NF_DEF_VAR (nid, "FOCE", NF_FLOAT, 2,dims, id_FOCE)
564        ierr = NF_PUT_ATT_TEXT (nid, id_FOCE, "title", 8,"Frac. Terre")
565        ierr = NF_DEF_VAR (nid, "FSIC", NF_FLOAT, 2,dims, id_FSIC)
566        ierr = NF_PUT_ATT_TEXT (nid, id_FSIC, "title", 8,"Frac. Terre")
567        ierr = NF_DEF_VAR (nid, "FLIC", NF_FLOAT, 2,dims, id_FLIC)
568        ierr = NF_PUT_ATT_TEXT (nid, id_FLIC, "title", 8,"Frac. Terre")
569!c
570        ierr = NF_ENDDEF(nid)
571!c
572
573! write the 'times'
574        do k=1,360
575#ifdef NC_DOUBLE
576          ierr = NF_PUT_VAR1_DOUBLE (nid,id_tim,k,DBLE(k))
577#else
578          ierr = NF_PUT_VAR1_REAL (nid,id_tim,k,FLOAT(k))
579#endif
580        enddo
581
582      endif ! of if (is_mpi_root.and.is_omp_root)
583
584! write the fields, after having collected them on master
585
586      call gather(phy_nat,phy_glo)
587      if (is_mpi_root.and.is_omp_root) then
588#ifdef NC_DOUBLE
589        ierr=NF_PUT_VAR_DOUBLE(nid,id_NAT,phy_glo)
590#else
591        ierr=NF_PUT_VAR_REAL(nid,id_NAT,phy_glo)
592#endif
593        if(ierr.ne.NF_NOERR) then
594          write(*,*) "writelim error with phy_nat"
595          write(*,*) NF_STRERROR(ierr)
596        endif
597      endif
598
599      call gather(phy_sst,phy_glo)
600      if (is_mpi_root.and.is_omp_root) then
601#ifdef NC_DOUBLE
602        ierr=NF_PUT_VAR_DOUBLE(nid,id_SST,phy_glo)
603#else
604        ierr=NF_PUT_VAR_REAL(nid,id_SST,phy_glo)
605#endif
606        if(ierr.ne.NF_NOERR) then
607          write(*,*) "writelim error with phy_sst"
608          write(*,*) NF_STRERROR(ierr)
609        endif
610      endif
611
612      call gather(phy_bil,phy_glo)
613      if (is_mpi_root.and.is_omp_root) then
614#ifdef NC_DOUBLE
615        ierr=NF_PUT_VAR_DOUBLE(nid,id_BILS,phy_glo)
616#else
617        ierr=NF_PUT_VAR_REAL(nid,id_BILS,phy_glo)
618#endif
619        if(ierr.ne.NF_NOERR) then
620          write(*,*) "writelim error with phy_bil"
621          write(*,*) NF_STRERROR(ierr)
622        endif
623      endif
624
625      call gather(phy_alb,phy_glo)
626      if (is_mpi_root.and.is_omp_root) then
627#ifdef NC_DOUBLE
628        ierr=NF_PUT_VAR_DOUBLE(nid,id_ALB,phy_glo)
629#else
630        ierr=NF_PUT_VAR_REAL(nid,id_ALB,phy_glo)
631#endif
632        if(ierr.ne.NF_NOERR) then
633          write(*,*) "writelim error with phy_alb"
634          write(*,*) NF_STRERROR(ierr)
635        endif
636      endif
637
638      call gather(phy_rug,phy_glo)
639      if (is_mpi_root.and.is_omp_root) then
640#ifdef NC_DOUBLE
641        ierr=NF_PUT_VAR_DOUBLE(nid,id_RUG,phy_glo)
642#else
643        ierr=NF_PUT_VAR_REAL(nid,id_RUG,phy_glo)
644#endif
645        if(ierr.ne.NF_NOERR) then
646          write(*,*) "writelim error with phy_rug"
647          write(*,*) NF_STRERROR(ierr)
648        endif
649      endif
650
651      call gather(phy_fter,phy_glo)
652      if (is_mpi_root.and.is_omp_root) then
653#ifdef NC_DOUBLE
654        ierr=NF_PUT_VAR_DOUBLE(nid,id_FTER,phy_glo)
655#else
656        ierr=NF_PUT_VAR_REAL(nid,id_FTER,phy_glo)
657#endif
658        if(ierr.ne.NF_NOERR) then
659          write(*,*) "writelim error with phy_fter"
660          write(*,*) NF_STRERROR(ierr)
661        endif
662      endif
663
664      call gather(phy_foce,phy_glo)
665      if (is_mpi_root.and.is_omp_root) then
666#ifdef NC_DOUBLE
667        ierr=NF_PUT_VAR_DOUBLE(nid,id_FOCE,phy_glo)
668#else
669        ierr=NF_PUT_VAR_REAL(nid,id_FOCE,phy_glo)
670#endif
671        if(ierr.ne.NF_NOERR) then
672          write(*,*) "writelim error with phy_foce"
673          write(*,*) NF_STRERROR(ierr)
674        endif
675      endif
676
677      call gather(phy_fsic,phy_glo)
678      if (is_mpi_root.and.is_omp_root) then
679#ifdef NC_DOUBLE
680        ierr=NF_PUT_VAR_DOUBLE(nid,id_FSIC,phy_glo)
681#else
682        ierr=NF_PUT_VAR_REAL(nid,id_FSIC,phy_glo)
683#endif
684        if(ierr.ne.NF_NOERR) then
685          write(*,*) "writelim error with phy_fsic"
686          write(*,*) NF_STRERROR(ierr)
687        endif
688      endif
689
690      call gather(phy_flic,phy_glo)
691      if (is_mpi_root.and.is_omp_root) then
692#ifdef NC_DOUBLE
693        ierr=NF_PUT_VAR_DOUBLE(nid,id_FLIC,phy_glo)
694#else
695        ierr=NF_PUT_VAR_REAL(nid,id_FLIC,phy_glo)
696#endif
697        if(ierr.ne.NF_NOERR) then
698          write(*,*) "writelim error with phy_flic"
699          write(*,*) NF_STRERROR(ierr)
700        endif
701      endif
702
703!  close file:
704      if (is_mpi_root.and.is_omp_root) then
705        ierr = NF_CLOSE(nid)
706      endif
707
708      end SUBROUTINE writelim
709
710!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
711
712      SUBROUTINE profil_sst(nlon, rlatd, type_profil, phy_sst)
713      use dimphy
714      IMPLICIT none
715!c
716      INTEGER nlon, type_profil, i, k, j
717      REAL :: rlatd(nlon), phy_sst(nlon, 360)
718      INTEGER imn, imx, amn, amx, kmn, kmx
719      INTEGER p, pplus, nlat_max
720      parameter (nlat_max=72)
721      REAL x_anom_sst(nlat_max)
722!c
723      if (klon.ne.nlon) stop'probleme de dimensions dans iniaqua'
724      do i=1,360
725!c      phy_sst(:,i) = 260.+50.*cos(rlatd(:))**2
726
727!c Rajout fbrlmd
728
729      if(type_profil.EQ.1)then
730!c     Méthode 1 "Control" faible plateau à l'Equateur
731      do j=1,klon
732       phy_sst(j,i)=273.+27.*(1-sin(1.5*rlatd(j))**2)
733!c        PI/3=1.047197551
734      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
735         phy_sst(j,i)=273.
736        endif
737      enddo
738      endif
739      if(type_profil.EQ.2)then
740!c     Méthode 2 "Flat" fort plateau à l'Equateur
741      do j=1,klon
742       phy_sst(j,i)=273.+27.*(1-sin(1.5*rlatd(j))**4)
743      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
744         phy_sst(j,i)=273.
745        endif
746      enddo
747      endif
748
749
750      if (type_profil.EQ.3) then
751!c     Méthode 3 "Qobs" plateau réel à l'Equateur
752      do j=1,klon
753        phy_sst(j,i)=273.+0.5*27.*(2-sin(1.5*rlatd(j))**2                            &
754       &                 -sin(1.5*rlatd(j))**4)
755      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
756         phy_sst(j,i)=273.
757        endif
758      enddo
759      endif
760
761      if (type_profil.EQ.4) then
762!c     Méthode 4 : Méthode 3 + SST+2 "Qobs" plateau réel à l'Equateur
763      do j=1,klon
764        phy_sst(j,i)=273.+0.5*29.*(2-sin(1.5*rlatd(j))**2                            &
765       &                 -sin(1.5*rlatd(j))**4)
766      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
767         phy_sst(j,i)=273.
768        endif
769      enddo
770      endif
771
772      if (type_profil.EQ.5) then
773!c     Méthode 5 : Méthode 3 + +2K "Qobs" plateau réel à l'Equateur
774      do j=1,klon
775        phy_sst(j,i)=273.+2.+0.5*27.*(2-sin(1.5*rlatd(j))**2                         &
776       &                 -sin(1.5*rlatd(j))**4)
777      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
778         phy_sst(j,i)=273.+2.
779        endif
780
781      enddo
782      endif
783
784      if(type_profil.EQ.6)then
785!c     Méthode 6 "cst" valeur constante de SST
786      do j=1,klon
787       phy_sst(j,i)=288.
788      enddo
789      endif
790
791
792        if(type_profil.EQ.7)then
793!c     Méthode 7 "cst" valeur constante de SST +2
794      do j=1,klon
795       phy_sst(j,i)=288.+2.
796      enddo
797      endif
798
799        p=0
800        if(type_profil.EQ.8)then
801!c     Méthode 8 profil anomalies SST du modÚle couplé AR4
802       do j=1,klon
803         if (rlatd(j).EQ.rlatd(j-1)) then
804       phy_sst(j,i)=273.+x_anom_sst(pplus)                                           &
805       &     +0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
806          if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
807            phy_sst(j,i)=273.+x_anom_sst(pplus)
808          endif
809        else
810          p=p+1
811          pplus=73-p
812        phy_sst(j,i)=273.+x_anom_sst(pplus)                                          &
813       &     +0.5*27.*(2-sin(1.5*rlatd(j))**2-sin(1.5*rlatd(j))**4)
814          if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
815            phy_sst(j,i)=273.+x_anom_sst(pplus)
816          endif
817          write (*,*) rlatd(j),x_anom_sst(pplus),phy_sst(j,i)
818        endif
819      enddo
820      endif
821
822      if (type_profil.EQ.9) then
823!c     Méthode 5 : Méthode 3 + -2K "Qobs" plateau réel à l'Equateur
824      do j=1,klon
825        phy_sst(j,i)=273.-2.+0.5*27.*(2-sin(1.5*rlatd(j))**2                         &
826       &                 -sin(1.5*rlatd(j))**4)
827      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
828         phy_sst(j,i)=273.-2.
829        endif
830      enddo
831      endif
832
833
834      if (type_profil.EQ.10) then
835!c     Méthode 10 : Méthode 3 + +4K "Qobs" plateau réel à l'Equateur
836      do j=1,klon
837        phy_sst(j,i)=273.+4.+0.5*27.*(2-sin(1.5*rlatd(j))**2                         &
838       &                 -sin(1.5*rlatd(j))**4)
839        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
840         phy_sst(j,i)=273.+4.
841        endif
842      enddo
843      endif
844
845      if (type_profil.EQ.11) then
846!c     Méthode 11 : Méthode 3 + 4CO2 "Qobs" plateau réel à l'Equateur
847      do j=1,klon
848        phy_sst(j,i)=273.+0.5*27.*(2-sin(1.5*rlatd(j))**2                            &
849       &                 -sin(1.5*rlatd(j))**4)
850        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
851         phy_sst(j,i)=273.
852        endif
853      enddo
854      endif
855
856      if (type_profil.EQ.12) then
857!c     Méthode 12 : Méthode 10 + 4CO2 "Qobs" plateau réel à l'Equateur
858      do j=1,klon
859        phy_sst(j,i)=273.+4.+0.5*27.*(2-sin(1.5*rlatd(j))**2                         &
860       &                 -sin(1.5*rlatd(j))**4)
861        if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
862         phy_sst(j,i)=273.+4.
863        endif
864      enddo
865      endif
866
867      if (type_profil.EQ.13) then
868!c     Méthode 13 "Qmax" plateau réel à l'Equateur augmenté !
869      do j=1,klon
870        phy_sst(j,i)=273.+0.5*29.*(2-sin(1.5*rlatd(j))**2                            &
871       &                 -sin(1.5*rlatd(j))**4)
872      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
873         phy_sst(j,i)=273.
874        endif
875      enddo
876      endif
877
878      if (type_profil.EQ.14) then
879!c     Méthode 13 "Qmax2K" plateau réel à l'Equateur augmenté +2K !
880      do j=1,klon
881        phy_sst(j,i)=273.+2.+0.5*29.*(2-sin(1.5*rlatd(j))**2                         &
882       &                 -sin(1.5*rlatd(j))**4)
883      if((rlatd(j).GT.1.0471975).OR.(rlatd(j).LT.-1.0471975))then
884         phy_sst(j,i)=273.
885        endif
886      enddo
887      endif
888
889      enddo
890
891!IM beg : verif profil SST: phy_sst
892       amn=MIN(phy_sst(1,1),1000.)
893       amx=MAX(phy_sst(1,1),-1000.)
894       DO k=1, 360
895       DO i=2, nlon
896        IF(phy_sst(i,k).LT.amn) THEN
897         amn=phy_sst(i,k)
898         imn=i
899         kmn=k
900        ENDIF
901        IF(phy_sst(i,k).GT.amx) THEN
902         amx=phy_sst(i,k)
903         imx=i
904         kmx=k
905        ENDIF
906       ENDDO
907       ENDDO
908!c
909       PRINT*,' debut avant writelim min max phy_sst',imn,kmn,amn,                   &
910       & imx,kmx,amx
911!IM end : verif profil SST: phy_sst
912
913       return
914       end SUBROUTINE profil_sst
915
916END MODULE phyaqua_mod
917
Note: See TracBrowser for help on using the repository browser.